C KTANGA SOURCE PV090527 26/04/30 21:15:47 12529 SUBROUTINE KTANGA(IPMOD0,IPCHE1,IPCHE2,IPCHE3,XPREC,DTPS,IKTSYM, & IPRIGI) *======================================================================= *= CALCUL DE LA MATRICE DE RIGIDITE TANGENTE = *======================================================================= *= Entrees : = *= --------- = *= IPMOD0 pointeur sur le mmodel = *= IPCHE1 pointeur sur le mchaml de contraintes = *= IPCHE2 pointeur sur le mchaml de variables internes = *= IPCHE3 pointeur sur le mchaml de caracteristiques = *= XPREC flottant precision = *= DTPS flottant pas de temps (modeles visco-plastiques) = *= IKTSYM =1 si matrice symetrique en sortie, =0 sinon = *======================================================================= *= Sortie : = *= -------- = *= IPRIGI pointeur sur matrice rigidite (=0 en cas d'erreur) = *======================================================================= *= Passage aux nouveaux chamelems par jm campenon le 05/91 = *= Mise a niveau FD/OF en 2009 = *======================================================================= IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC CCGEOME -INC SMCHAML -INC SMCOORD -INC SMELEME -INC SMINTE -INC SMLREEL -INC SMMODEL -INC SMRIGID INTEGER OOOVAL -INC TMPTVAL SEGMENT WRK1 REAL*8 DDHOOK(NSTRS,NSTRS),DDHOMU(NSTRS,NSTRS), & REL(LRE,LRE),XE(3,NBBB) ENDSEGMENT SEGMENT WRK2 REAL*8 SHPWRK(6,NBNO),BGENE(NSTRS,LRE) ENDSEGMENT SEGMENT WRK3 REAL*8 WORK(LW) ENDSEGMENT SEGMENT WRK4 REAL*8 BPSS(3,3),XEL(3,NBBB) ENDSEGMENT SEGMENT WRK5 INTEGER NTRAC1,NTRAC2 ENDSEGMENT * POUR LES MATERIAUX a "TROPIE" (PASSAGE DE LA MATRICE DE ROTATION) SEGMENT WTRAV REAL*8 TXR(IDIM,IDIM) ENDSEGMENT SEGMENT NOTYPE CHARACTER*16 TYPE(NBTYPE) ENDSEGMENT C- Nombre de points maximal pour stocker une courbe de traction PARAMETER (LTRAC=2*75) * INTTYP correspond au type de points d'integration utilise par KTAN PARAMETER ( INTTYP=3 ) DIMENSION TRAC(LTRAC) DIMENSION CRIGI(12),CMASS(12) DIMENSION A(4,60),BB(3,60),PP(4,4) * Petit tableau des "couleurs" des relations de conformite DIMENSION LCOLOR(6) DATA LCOLOR / 1, 3, 6, 10, 16, 24 / CHARACTER*8 CMATE CHARACTER*(NCONCH) CONM PARAMETER ( NINF=3 ) INTEGER INFOS(NINF) LOGICAL lsupma, BDPGE,BPLAN,BMATE *======================================================================* *= 1 - INITIALISATIONS ET VERIFICATIONS =* *======================================================================* bmate =.FALSE. IPRIGI=0 KERRE=0 * Modele "deroule" (uniquement "MECANIQUE", "LIQUIDE" ou "POREUX") CALL PIMODL(IPMOD0,IPMODL,MAILDG,0) IF (IPMODL.EQ.0) RETURN * Verification du support du mchaml de CONTRAINTES CALL REDUAF(IPCHE1,IPMOD0,IPCH_Z,0,IRET,KERRE) IF (IRET.NE.1) CALL ERREUR(KERRE) IF (IERR.NE.0) GOTO 550 IPCHE1=IPCH_Z CALL QUESUP(IPMODL,IPCHE1,INTTYP,0,ISUPCO,IRET) IF (ISUPCO.GT.1) GOTO 550 * Verification du support du mchaml de VARIABLES INTERNES CALL REDUAF(IPCHE2,IPMOD0,IPCH_Z,0,IRET,KERRE) IF (IRET.NE.1) CALL ERREUR(KERRE) IF (IERR.NE.0) GOTO 550 IPCHE2=IPCH_Z CALL QUESUP(IPMODL,IPCHE2,INTTYP,0,ISUPVA,IRET) IF (ISUPVA.GT.1) GOTO 550 * Verification du support du mchaml de CARACTERISTIQUES CALL REDUAF(IPCHE3,IPMOD0,IPCH_Z,0,IRET,KERRE) IF (IRET.NE.1) CALL ERREUR(KERRE) IF (IERR.NE.0) GOTO 550 IPCHE3=IPCH_Z CALL QUESUP(IPMODL,IPCHE3,INTTYP,0,ISUPMA,IRET) IF (ISUPMA.GT.1) GOTO 550 * Activation du modele MMODEL=IPMODL NSOUS=KMODEL(/1) * Initialisations de l'objet RIGIDITE "matrice tangente" NRIGEL=NSOUS SEGINI,MRIGID IPRIGI=MRIGID MTYMAT='RIGIDITE' ICHOLE=0 IFORIG=IFOUR IMGEO1=0 IMGEO2=0 ISUPEQ=0 NHRM=NIFOUR melpha = 0 * Indicateur de mode de calcul en 2D plan BPLAN = IFOUR.EQ.-2 .OR. IFOUR.EQ.-1 .OR. IFOUR.EQ.-3 * Type des composantes NBTYPE=1 SEGINI,NOTYPE TYPE(1)='REAL*8' MOTYR8=NOTYPE *======================================================================* *= 2 - BOUCLE SUR LES SOUS-ZONES DU MODELE (Fin = etiquette 500) =* *======================================================================* ISOU = 0 DO 500 ISOUS=1,NSOUS IMODEL=KMODEL(ISOUS) IPMOD1=IMODEL *----------------------------------------------------------------------- *- 2.1 - Intialisations et activations de segments *----------------------------------------------------------------------- MELE = imodel.NEFMOD IPMAIL = imodel.IMAMOD IIPDPG = imodel.IPDPGE IIPDPG = IPTPOI(IIPDPG) IPINF = 0 * Cas particulier des relations de conformites IF (MELE.EQ.22) GOTO 5001 IF (MELE.EQ.259) GOTO 5001 * Verifications sur la formulation CONM = CONMOD CMATE = CMATEE MATE = IMATEE MAPL = INATUU BMATE = (CMATE.EQ.'UNIDIREC').OR.(CMATE.EQ.'ORTHOTRO').OR. & (CMATE.EQ.'ANISOTRO') * Information sur l'element fini IF (INFMOD(/1).LT.2+INTTYP) THEN write(ioimp,*) 'KTANGA - INFMOD(/1) =',infmod(/1),'<',2+inttyp call erreur(5) ENDIF NBGS =INFELE(4) C* ICARA=INFELE(5) NBPGAU=INFELE(6) LW =INFELE(7) IPORE=INFELE(8) LRE =INFELE(9) LHOOK=INFELE(10) MFR =INFELE(13) IELE =INFELE(14) NDDL =INFELE(15) C* NSTRS=INFELE(16) IPMINT=INFMOD(2+INTTYP) IPMIN1=INFMOD(3) MINTE=IPMINT IF (MFR.EQ.57.OR.MFR.EQ.59) THEN IPPORE = NBNNE(NUMGEO(MELE)) IF (IFOUR.EQ.1.OR.IFOUR.EQ.-3) THEN LHOOK=6 ELSE LHOOK=4 ENDIF ELSE IF (MFR.EQ.33) THEN IPPORE = NBNNE(NUMGEO(MELE)) ELSE IPPORE = 0 ENDIF C Coque integree ou non ? NPINT = INFMOD(1) CALL INFDPG(MFR,IFOUR, BDPGE,NDPGE) * Coordonnees du point support des deformations planes generalisees IF (BDPGE) THEN IF (IIPDPG.LE.0) THEN CALL ERREUR(925) CALL ERREUR(5) GOTO 551 ENDIF IREF=(IIPDPG-1)*(IDIM+1) XDPGE=XCOOR(IREF+1) YDPGE=XCOOR(IREF+2) ELSE XDPGE = 0.D0 YDPGE = 0.D0 ENDIF *----------------------------------------------------------------------- *- 2.2 - Preparation des objets resultats DESCR et XMATRI *----------------------------------------------------------------------- * Si necessaire PARTITIONNEMENT du segment XMATRI 5001 CONTINUE LTRK=OOOVAL(1,4) IF (LTRK.EQ.0) LTRK=OOOVAL(1,1) LTRK=MAX(LTRK,2**24) IPT1=IPMAIL SEGACT,IPT1 NBNN1 =IPT1.NUM(/1) NBELE1=IPT1.NUM(/2) IF (MELE.EQ.22) LRE=NBNN1 IF (MELE.EQ.259) LRE=NBNN1 * Traitements particuliers pour penalisation milieu poreux IDECAP = 0 IF (MELE.GE.79.AND.MELE.LE.83) THEN IDECAP = 1 LRE = LRE + 2*NBNN1 - IPORE ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN IDECAP=1 LRE = LRE + (3*NBNN1 - IPORE)/2 - NBSOM(IELE) ELSE IF (MELE.GE.173.AND.MELE.LE.177) THEN IDECAP=2 LRE = LRE + (2*NBNN1 - IPORE)*IDECAP ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN IDECAP=3 LRE = LRE + (2*NBNN1 - IPORE)*IDECAP ENDIF * Ajout a la taille en mots de la matrice des infos du segment LSEG=LRE*LRE*NBELE1 + 16 NBLPRT=(LSEG-1)/LTRK+1 NBLMAX=(NBELE1-1)/NBLPRT+1 NBLPRT=(NBELE1-1)/NBLMAX+1 * write(ioimp,*) ' ktanga nblprt nblmax ',NBLPRT,NBLMAX,NBELE1 MELEME=IPT1 * Boucle (5000) de PARTITIONNEMENT du segment XMATRI DO 5000 IPRT = 1,NBLPRT ISOU=ISOU+1 IF (ISOU.GT.IRIGEL(/2)) THEN NRIGEL=ISOU SEGADJ,MRIGID ENDIF IF (NBLPRT.GT.1) THEN JPRT=(IPRT-1)*NBLMAX SEGACT,IPT1 NBSOUS=0 NBREF=0 NBNN=NBNN1 NBELEM=MIN(NBLMAX,NBELE1-JPRT) * write(ioimp,*) ' creation segment ',nbnn,nbelem SEGINI,MELEME ITYPEL=IPT1.ITYPEL DO I=1,NBELEM IB=I+JPRT DO J=1,NBNN NUM(J,I)=IPT1.NUM(J,IB) ENDDO ICOLOR(I)=IPT1.ICOLOR(I) ENDDO ENDIF IPMAIL=MELEME * Fin du traitement particulier en cas de PARTITIONNEMENT du XMATRI * Quelques initialisations suite au partionnement IPDSCR = 0 IPMADG = 0 IPMATR = 0 IRIGE7 = 0 NMATR = 0 NMATF = 0 IVAMAT = 0 NCARA = 0 NCARF = 0 IVACAR = 0 NVARI = 0 NVARF = 0 IVARI = 0 IVACON = 0 * Activation du MELEME support des rigidites MELEME=IPMAIL NBNN =NUM(/1) NBELEM=NUM(/2) * Cas particulier des relations de conformites IF (MELE.EQ.22) GOTO 22 IF (MELE.EQ.259) GOTO 259 * Modification du MELEME pour les deformations planes generalisees IF (BDPGE) THEN NBNA=NBNN NBNN=NBNA+1 NBREF=0 NBSOUS=0 SEGINI,IPT2 IPT2.ITYPEL=28 DO I=1,NBELEM DO J=1,NBNA IPT2.NUM(J,I)=NUM(J,I) ENDDO IPT2.NUM(NBNN,I)=IIPDPG IPT2.ICOLOR(I)=ICOLOR(I) ENDDO IPMAGD=IPT2 ENDIF * Recherche des noms d'inconnues primales et duales NOMID=LNOMID(1) if (nomid.eq.0) then write(ioimp,*) 'KTANGA : MODEPL = lnomid(1) = 0' call erreur(5) return endif MODEPL = NOMID NDEPL = nomid.LESOBL(/2) c* nfac = nomid.LESFAC(/2) NOMID=LNOMID(2) if (nomid.eq.0) then write(ioimp,*) 'KTANGA : MOFORC = lnomid(2) = 0' call erreur(5) return endif MOFORC = NOMID NFORC = nomid.LESOBL(/2) c* nfac = nomid.LESFAC(/2) IF (NDEPL.EQ.0.OR.NFORC.EQ.0.OR.NDEPL.NE.NFORC) THEN call erreur(5) return ENDIF * Initialisation du segment DESCR NLIGRP = LRE NLIGRD = LRE SEGINI,DESCR IPDSCR=DESCR * Remplissage du segment DESCRipteur NCOMP = NDEPL NBNNS = NBNN IF (MFR.EQ.33.OR.MFR.EQ.57.OR.MFR.EQ.59) THEN NCOMP = NDEPL-IDECAP ENDIF IF (BDPGE) THEN NCOMP = NDEPL-NDPGE NBNNS = NBNN-1 ENDIF IF (MFR.EQ.19.OR.MFR.EQ.21) NBNNS=NBNN/2 IDDL=1 * Cas du macro-element IF (MFR.EQ.61)THEN DO i=1,3 NOELEP(i )=1 NOELEP(i+3)=3 ENDDO NOELEP(7)=2 NOELEP(8)=2 DO i=1,LRE NOELED(i)=NOELEP(i) ENDDO NOMID=MODEPL DO i=1,3 LISINC(i )=LESOBL(i) LISINC(i+3)=LESOBL(i) ENDDO LISINC(7)=LESOBL(4) LISINC(8)=LESOBL(5) NOMID=MOFORC DO i=1,3 LISDUA(i )=LESOBL(i) LISDUA(i+3)=LESOBL(i) ENDDO LISDUA(7)=LESOBL(4) LISDUA(8)=LESOBL(5) * Cas general ELSE * Erreur dans les dimensions de DESCR (mode de calcul incorrect) IF (NBNNS*NCOMP.GT.NLIGRD) THEN KERRE=717 GOTO 515 ENDIF NDUM=NBNNS IF (MELE.GE.108.AND.MELE.LE.110) THEN NFAC=(3*NBNN-IPORE)/2 NDUM=MIN(NBNNS,NFAC) ENDIF DO INOEUD=1,NDUM DO ICOMP=1,NCOMP NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODEPL LISINC(IDDL)=LESOBL(ICOMP) NOMID=MOFORC LISDUA(IDDL)=LESOBL(ICOMP) IDDL=IDDL+1 ENDDO ENDDO ENDIF * Cas particulier des deformations planes generalisees IF (BDPGE) THEN DO ICOMP=(NDPGE-1),0,-1 NOELEP(IDDL)=NBNN NOELED(IDDL)=NBNN NOMID=MODEPL LISINC(IDDL)=LESOBL(NDEPL-ICOMP) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NFORC-ICOMP) IDDL=IDDL+1 ENDDO ENDIF * Cas particulier des milieux poreux (pression aux sommets en 1er) IF (MFR.EQ.33) THEN DO INOEUD=1,NBSOM(IELE) NOELEP(IDDL)=IBSOM(NSPOS(IELE)+INOEUD-1) NOELED(IDDL)=NOELEP(IDDL) NOMID=MODEPL LISINC(IDDL)=LESOBL(NDEPL) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDEPL) IDDL=IDDL+1 ENDDO IF (MELE.GE.79.AND.MELE.LE.83) THEN DO 1105 INOEUD=1,NBNN DO i=1,NBSOM(IELE) IF (INOEUD.EQ.IBSOM(NSPOS(IELE)+i-1)) GOTO 1105 ENDDO NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODEPL LISINC(IDDL)=LESOBL(NDEPL) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDEPL) IDDL=IDDL+1 1105 CONTINUE ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN DO INOEUD=NFAC+1,NBNN NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODEPL LISINC(IDDL)=LESOBL(NDEPL) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDEPL) IDDL=IDDL+1 ENDDO DO 1110 INOEUD=1,NFAC DO i=1,NBSOM(IELE) IF(INOEUD.EQ.IBSOM(NSPOS(IELE)+i-1)) GOTO 1110 ENDDO NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODEPL LISINC(IDDL)=LESOBL(NDEPL) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDEPL) IDDL=IDDL+1 1110 CONTINUE ENDIF ELSE IF (MFR.EQ.57.OR.MFR.EQ.59) THEN IF (MELE.GE.173.AND.MELE.LE.182) THEN DO IPR=1,IDECAP NDECAP = NDEPL-IDECAP+IPR DO INOEUD=1,NBSOM(IELE) NOELEP(IDDL)=IBSOM(NSPOS(IELE)+INOEUD-1) NOELED(IDDL)=NOELEP(IDDL) NOMID=MODEPL LISINC(IDDL)=LESOBL(NDECAP) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDECAP) IDDL=IDDL+1 ENDDO DO 1205 INOEUD=1,NBNN DO i=1,NBSOM(IELE) IF (INOEUD.EQ.IBSOM(NSPOS(IELE)+i-1)) GOTO 1205 ENDDO NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODEPL LISINC(IDDL)=LESOBL(NDECAP) NOMID=MOFORC LISDUA(IDDL)=LESOBL(NDECAP) IDDL=IDDL+1 1205 CONTINUE ENDDO ENDIF * Cas des elements raccord ELSE IF (MFR.EQ.19.OR.MFR.EQ.21) THEN CALL IDPRIM(IMODEL,MFR+1000,MODPL,NDEPL,ndum) NOMID=MODPL SEGACT,NOMID CALL IDDUAL(IMODEL,MFR+1000,MOFRC,NFORC,ndum) NOMID=MOFRC SEGACT,NOMID DO INOEUD=NBNNS+1,NBNN DO ICOMP=1,NDEPL NOELEP(IDDL)=INOEUD NOELED(IDDL)=INOEUD NOMID=MODPL LISINC(IDDL)=LESOBL(ICOMP) NOMID=MOFRC LISDUA(IDDL)=LESOBL(ICOMP) IDDL=IDDL+1 ENDDO ENDDO NOMID=MODPL SEGSUP,NOMID NOMID=MOFRC SEGSUP,NOMID ENDIF * Initialisation du segment XMATRI contenant les matrices elementaires * de la sous-zone (NBELEM = nombre d'elements dans la sous-zone =MELEME) NELRIG=NBELEM RIGREL=0 SEGINI,XMATRI IPMATR=XMATRI * Quelques donnes utiles pour le segment MRIGID IF (BDPGE) THEN ** MELEME=IPMAIL <- MELEME segment actif et pointe tjs sur IPMAIL NBNN=NUM(/1) ELSE IPMAGD=IPMAIL ENDIF IF (MAPL.EQ.35.OR.MAPL.EQ.54.OR.MAPL.EQ.56.OR.MAPL.EQ.111) THEN IRIGE7=2 ELSE IF (MFR.EQ.57.OR.MFR.EQ.59) THEN IRIGE7=2 ELSE IRIGE7=0 ENDIF * En cas de rendement IRIGE7=2 (cf. RIGI1.ESO) *----------------------------------------------------------------------- *- 2.3 - Analyse des champs par element fournis en entree *----------------------------------------------------------------------- * Creation du tableau infos (contraintes, variables internes) CALL IDENT(IPMAIL,CONM,IPCHE2,IPCHE1,INFOS,IRET) IF (IRET.EQ.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF * Recherche des noms de composantes du champ de CONTRAINTEs NOMID=LNOMID(4) if (nomid.eq.0) then write(ioimp,*) 'KTANGA : MOCONT = lnomid(4) = 0' call erreur(5) return endif MOCONT = NOMID NSTRS = nomid.LESOBL(/2) C* nfac = nomid.LESFAC(/2) * Verification de leur presence CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOCONT,MOTYR8,1,INFOS,3,IVACON) IF (IERR.NE.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF IF (ISUPCO.EQ.1) THEN CALL VALCHE(IVACON,NSTRS,IPMINT,IPPORE,MOCONT,MELE) IF (IERR.NE.0) THEN ISUPCO=0 KERRE=-ABS(IERR) GOTO 515 ENDIF ENDIF * Recherche des noms de composantes du champ des variables internes NOMID=LNOMID(10) if (nomid.eq.0) then write(ioimp,*) 'KTANGA : MOVARI = lnomid(10) = 0' KERRE=76 MOTERR(1:4)='VARI' MOTERR(5:8)=NOMTP(MELE) GOTO 515 endif MOVARI = NOMID NVARI = nomid.LESOBL(/2) NVARF = nomid.LESFAC(/2) NVART=NVARI+NVARF * Type des composantes notype = motyr8 IF (CMATE.EQ.'SECTION') THEN NBTYPE=1 SEGINI,NOTYPE TYPE(1)='POINTEURMCHAML ' ENDIF MOTYPE=NOTYPE * Verification de leur presence CALL KOMCHA(IPCHE2,IPMAIL,CONM,MOVARI,MOTYPE,1,INFOS,3,IVARI) IF (MOTYPE.NE.MOTYR8) SEGSUP,NOTYPE IF (IERR.NE.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF IF (ISUPVA.EQ.1) THEN CALL VALCHE(IVARI,NVART,IPMINT,IPPORE,MOVARI,MELE) IF (IERR.NE.0) THEN ISUPVA=0 KERRE=-ABS(IERR) GOTO 515 ENDIF ENDIF * Creation du tableau infos (variables internes, caracteristiques) CALL IDENT(IPMAIL,CONM,IPCHE2,IPCHE3,INFOS,IRET) IF (IRET.EQ.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF * Recuperation des noms de composantes de caracteristiques materielles * Sauf cas particulier TYPE='REAL*8' NBROBL=0 NBRFAC=0 NOMID=0 lsupma=.TRUE. NOTYPE=MOTYR8 * Element de barre et Acier Unidirirectionnel IF (MAPL.EQ.40.AND.MFR.EQ.27) THEN NBROBL=1 SEGINI,NOMID LESOBL(1)='YOUN' ELSE IF (MFR.EQ.7.AND.CMATE.EQ.'SECTION') THEN NBROBL=2 NBRFAC=1 SEGINI,NOMID LESOBL(1)='MODS' LESOBL(2)='MATS' LESFAC(1)='MAHO' NBTYPE=3 SEGINI,NOTYPE TYPE(1)='POINTEURMMODEL ' TYPE(2)='POINTEURMCHAML ' TYPE(3)='POINTEURLISTREEL' * Cas POI1 -- MODAL -- MFR=26 ==> pas traite dans la suite ELSE IF (MFR.EQ.26) THEN NBROBL=3 SEGINI,NOMID LESOBL(1)='FREQ' LESOBL(2)='MASS' LESOBL(3)='DEFO' NBTYPE=3 SEGINI,NOTYPE TYPE(1)='REAL*8' TYPE(2)='REAL*8' TYPE(3)='POINTEURCHPOINT' * Cas POI1 -- STATIQUE -- MFR=28 ==> pas traite dans la suite ELSE IF (MFR.EQ.28) THEN NBROBL=3 SEGINI,NOMID LESOBL(1)='DEFO' LESOBL(2)='RIDE' LESOBL(3)='MADE' NBTYPE=1 SEGINI,NOTYPE TYPE(1)='POINTEURCHPOINT' * Cas Orthotrope, Anisotrope et Unidirectionnel ELSE IF (BMATE) THEN * Materiau Unidirirectionnel C*? IF (FORMOD(/1).EQ.'MECANIQUE'.AND.CMATE.EQ.'UNIDIREC') THEN IF (CMATE.EQ.'UNIDIREC') THEN IF (MFR.EQ.1.AND.IDIM.EQ.3) THEN NBROBL=7 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='V1X ' LESOBL(3)='V1Y ' LESOBL(4)='V1Z ' LESOBL(5)='V2X ' LESOBL(6)='V2Y ' LESOBL(7)='V2Z ' ELSE NBROBL=3 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='V1X ' LESOBL(3)='V1Y ' ENDIF * Materiau orthotrope plastique 'ECROUIS_DECOU' ELSE IF (CMATE.EQ.'ORTHOTRO'.AND.MAPL.EQ.67) THEN NBROBL=6 SEGINI,NOMID LESOBL(1)='YG1 ' LESOBL(2)='YG2 ' LESOBL(3)='NU12' LESOBL(4)='G12 ' LESOBL(5)='V1X ' LESOBL(6)='V1Y ' * Autres Materiaux orthotropes et anisotropes ELSE IF (LNOMID(6).NE.0) THEN lsupma=.FALSE. NOMID=LNOMID(6) NBROBL=LESOBL(/2) NBRFAC=LESFAC(/2) ELSE CALL IDMATR(MFR,IPMOD1,MOMATR,NBROBL,NBRFAC) ENDIF * Cas particulier : Mistral (10 composantes = listes de reels) IF (MAPL.EQ.94) THEN NBTYPE=NBROBL+NBRFAC SEGINI,NOTYPE DO i=1,NBTYPE TYPE(i)='REAL*8' ENDDO NLDEB=NBROBL-9 DO i=NLDEB,NBROBL TYPE(i)='POINTEURLISTREEL' ENDDO ENDIF ENDIF * Materiaux ISOTROPEs ELSE IF (CMATE.EQ.'ISOTROPE') THEN IF (MFR.EQ.35) THEN IF (MAPL.EQ.35) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='KS ' LESOBL(2)='KN ' LESOBL(3)='PHI ' LESOBL(4)='MU ' ELSE NBROBL=2 SEGINI,NOMID LESOBL(1)='KS ' LESOBL(2)='KN ' ENDIF * Element joint cisaillement 2D ELSE IF (MFR.EQ.53) THEN NBROBL=1 SEGINI,NOMID LESOBL(1)='KS ' * Elements POREUX isotropes ELSE IF (FORMOD(1).EQ.'POREUX') THEN IF (MELE.GE.79.AND.MELE.LE.83) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='COB ' LESOBL(4)='MOB ' ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='KS ' LESOBL(2)='KN ' LESOBL(3)='COB ' LESOBL(4)='MOB ' ELSE IF (MELE.GE.173.AND.MELE.LE.177) THEN NBROBL=10 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='COP1' LESOBL(4)='COP2' LESOBL(5)='CPP1' LESOBL(6)='CPP2' LESOBL(7)='KK11' LESOBL(8)='KK12' LESOBL(9)='KK21' LESOBL(10)='KK22' ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN NBROBL=17 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='COP1' LESOBL(4)='COP2' LESOBL(5)='COP3' LESOBL(6)='CPP1' LESOBL(7)='CPP2' LESOBL(8)='CPP3' LESOBL(9)='KK11' LESOBL(10)='KK12' LESOBL(11)='KK13' LESOBL(12)='KK21' LESOBL(13)='KK22' LESOBL(14)='KK23' LESOBL(15)='KK31' LESOBL(16)='KK32' LESOBL(17)='KK33' ENDIF ELSE IF (MAPL.EQ.1) THEN NBROBL=3 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='SIGY' ELSE IF (MAPL.EQ.2.OR.MAPL.EQ.14) THEN NBROBL=2 NBRFAC=2 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESFAC(1)='SIGF' LESFAC(2)='TRAC' NBTYPE=4 SEGINI,NOTYPE TYPE(1)='REAL*8' TYPE(2)='REAL*8' TYPE(3)='REAL*8' TYPE(4)='POINTEUREVOLUTIO' ELSE IF (MAPL.EQ.3) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='LTR ' LESOBL(4)='LCS ' ELSE IF (MAPL.EQ.4) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='SIGY' LESOBL(4)='H ' ELSE IF (MAPL.EQ.5) THEN NBROBL=3 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='TRAC' NBTYPE=3 SEGINI,NOTYPE TYPE(1)='REAL*8' TYPE(2)='REAL*8' TYPE(3)='POINTEUREVOLUTIO' * Modele Drucker Prager ELSE IF (MAPL.EQ.15) THEN NBROBL=11 SEGINI,NOMID LESOBL(1) ='YOUN' LESOBL(2) ='NU ' LESOBL(3) ='ETA ' LESOBL(4) ='MU ' LESOBL(5) ='KL ' LESOBL(6) ='GAMM' LESOBL(7) ='DELT' LESOBL(8) ='ALFA' LESOBL(9) ='BETA' LESOBL(10)='K ' LESOBL(11)='H ' * Modele visco-plastique parfait ELSE IF (MAPL.EQ.43) THEN NBROBL=5 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='SIGY' LESOBL(4)='N ' LESOBL(5)='K ' * Modele Betocyclique ELSE IF (MAPL.EQ.54) THEN NBROBL=13 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='HHH1' LESOBL(4)='FTPE' LESOBL(5)='FCPE' LESOBL(6)='FTGR' LESOBL(7)='FCGR' LESOBL(8)='EPSO' LESOBL(9)='WOR0' LESOBL(10)='LCAT' LESOBL(11)='LCAC' LESOBL(12)='TREV' LESOBL(13)='COEV' NBTYPE=13 SEGINI,NOTYPE DO i=1,NBTYPE-2 TYPE(i)='REAL*8' ENDDO TYPE(12)='POINTEUREVOLUTIO' TYPE(13)='POINTEUREVOLUTIO' * Rotating Crack ELSE IF (MAPL.EQ.55) THEN NBROBL=6 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='FTRA' LESOBL(4)='EPSR' LESOBL(5)='FRES' LESOBL(6)='BETA' * BCN-MRS-Lade (MAPL=111) ELSE IF (MAPL.EQ.111) THEN NBROBL=20 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='PC ' LESOBL(4)='PA ' LESOBL(5)='QA ' LESOBL(6)='EXPM' LESOBL(7)='E ' LESOBL(8)='K1 ' LESOBL(9)='K2 ' LESOBL(10)='ETAB' LESOBL(11)='EXPV' LESOBL(12)='EPSI' LESOBL(13)='N ' LESOBL(14)='CCON' LESOBL(15)='EXPL' LESOBL(16)='PCAP' LESOBL(17)='EXPR' LESOBL(18)='CCAP' LESOBL(19)='PHI ' LESOBL(20)='ALP ' * BCN-J2 (MAPL=112) ELSE IF (MAPL.EQ.112) THEN NBROBL=6 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='SIG0' LESOBL(4)='SIGI' LESOBL(5)='KISO' LESOBL(6)='VELO' * BCN-Rounded Hyperbolic Mohr-Coulomb (MAPL=113) ELSE IF (MAPL.EQ.113) THEN NBROBL=4 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' LESOBL(3)='COHE' LESOBL(4)='PHI ' * Autres modeles ISOTROPEs : elasticite ELSE NBROBL=2 SEGINI,NOMID LESOBL(1)='YOUN' LESOBL(2)='NU ' ENDIF * Autres cas ? ELSE IF (LNOMID(6).NE.0) THEN lsupma=.FALSE. NOMID=LNOMID(6) NBROBL=LESOBL(/2) NBRFAC=LESFAC(/2) ELSE CALL IDMATR(MFR,IPMOD1,MOMATR,NBROBL,NBRFAC) ENDIF IF (CMATE.EQ.'SECTION') THEN NBTYPE=3 SEGINI,NOTYPE TYPE(1)='POINTEURMMODEL' TYPE(2)='POINTEURMCHAML' TYPE(3)='POINTEURLISTREEL' ENDIF ENDIF MOMATR=NOMID MOTYPE=NOTYPE IF (MOMATR.EQ.0) THEN if (motype.NE.MOTYR8) SEGSUP,NOTYPE KERRE=591 GOTO 515 ENDIF CALL KOMCHA(IPCHE3,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) if (motype.NE.MOTYR8) SEGSUP,NOTYPE IF (lsupma) SEGSUP,NOMID IF (IERR.NE.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF NMATR=NBROBL NMATF=NBRFAC NMATT=NMATR+NMATF IF (ISUPMA.EQ.1) THEN CALL VALCHE(IVAMAT,NMATT,IPMINT,IPPORE,MOMATR,MELE) IF (IERR.NE.0) THEN ISUPMA=0 KERRE=-ABS(IERR) GOTO 515 ENDIF ENDIF * Recuperation des noms de composantes de caracteristiques geometriques * Sauf cas particulier MOTYPE = segment NBTYPE=1 et TYPE(1)='REAL*8' NOTYPE=MOTYR8 NBROBL=0 NBRFAC=0 NOMID=0 IVECT=0 * Massif ou certains elements poreux en contraintes planes IF ( (MFR.EQ.1 .OR. MFR.EQ.31 .OR. & (MELE.GE.79.AND.MELE.LE.83) .OR. & (MELE.GE.173.AND.MELE.LE.182) ) & .AND. IFOUR.EQ.-2 ) THEN NBRFAC=1 SEGINI,NOMID LESFAC(1)='DIM3' * Cas des coques ELSE IF (MFR.EQ.3 .OR. MFR.EQ.5 .OR. MFR.EQ.9) THEN NBROBL=1 IF (MFR.EQ.3.AND.IFOUR.EQ.-2) THEN NBRFAC=3 ELSE NBRFAC=2 ENDIF SEGINI,NOMID LESOBL(1)='EPAI' LESFAC(1)='EXCE' IF (MFR.EQ.3.AND.IFOUR.EQ.-2) LESFAC(2)='DIM3' LESFAC(NBRFAC)='CALF' * Donnees pour les poutres ELSE IF (MFR.EQ.7) THEN IF (CMATE.NE.'SECTION' ) THEN IF (BPLAN) THEN NBRFAC=1 NBROBL=2 SEGINI,NOMID LESOBL(1)='SECT' LESOBL(2)='INRZ' LESFAC(1)='SECY' ELSE NBROBL=4 NBRFAC=5 SEGINI,NOMID LESOBL(1)='TORS' LESOBL(2)='INRY' LESOBL(3)='INRZ' LESOBL(4)='SECT' LESFAC(1)='SECY' LESFAC(2)='SECZ' LESFAC(3)='VX' LESFAC(4)='VY' LESFAC(5)='VZ' IVECT=1 ENDIF ELSE NBRFAC=3 SEGINI,NOMID LESFAC(1)='VX' LESFAC(2)='VY' LESFAC(3)='VZ' IVECT=1 ENDIF C Donnees pour les TUYAUX ELSE IF (MFR.EQ.13) THEN NBROBL=2 NBRFAC=6 SEGINI,NOMID LESOBL(1)='EPAI' LESOBL(2)='RAYO' LESFAC(1)='RACO' LESFAC(2)='PRES' LESFAC(3)='CISA' LESFAC(4)='VX' LESFAC(5)='VY' LESFAC(6)='VZ' IVECT=1 C Donnees pour le LINESPRING ELSE IF (MFR.EQ.15) THEN NBROBL=5 SEGINI,NOMID LESOBL(1)='EPAI' LESOBL(2)='FISS' LESOBL(3)='VX ' LESOBL(4)='VY ' LESOBL(5)='VZ ' C Donnees pour le TUYAU FISSURE ELSE IF (MFR.EQ.17) THEN NBROBL=9 SEGINI,NOMID LESOBL(1)='RAYO' LESOBL(2)='EPAI' LESOBL(3)='VX ' LESOBL(4)='VY ' LESOBL(5)='VZ ' LESOBL(6)='VXF ' LESOBL(7)='VYF ' LESOBL(8)='VZF ' LESOBL(9)='ANGL' * Section pour les barres - uniaxial ELSE IF (MFR.EQ.27 .AND. CMATE.NE.'NODAL') THEN NBROBL=1 SEGINI,NOMID LESOBL(1)='SECT' * Elements homogeneises ELSE IF (MFR.EQ.37) THEN IF (IFOUR.EQ.1.OR.IFOUR.EQ.0.OR.IFOUR.EQ.2) THEN NBROBL=5 SEGINI,NOMID LESOBL(1)='SCEL' LESOBL(2)='SFLU' LESOBL(3)='EPS ' LESOBL(4)='SECT' LESOBL(5)='INRZ ' ELSE NBROBL=3 SEGINI,NOMID LESOBL(1)='SCEL' LESOBL(2)='SFLU' LESOBL(3)='EPS ' ENDIF * Element TUYO ELSE IF (MFR.EQ.39) THEN NBROBL=2 NBRFAC=5 SEGINI,NOMID LESOBL(1)='EPAI' LESOBL(2)='RAYO' LESFAC(1)='RACO' LESFAC(2)='PRES' LESFAC(3)='VX' LESFAC(4)='VY' LESFAC(5)='VZ' IVECT=1 * Element tuyau acoustique pure ELSE IF (MFR.EQ.41) THEN NBROBL=1 NBRFAC=1 SEGINI,NOMID LESOBL(1)='RAYO' LESFAC(1)='RACO' * Donnees pour les barres excentrees ELSE IF (MFR.EQ.49) THEN NBROBL=6 SEGINI,NOMID LESOBL(1)='SECT' LESOBL(2)='EXCZ' LESOBL(3)='EXCY' LESOBL(4)='VX ' LESOBL(5)='VY ' LESOBL(6)='VZ ' * Donnees geometriques pour l'element LIA2 de liaison a 2 noeuds ELSE IF (MFR.EQ.51) THEN NBROBL=9 SEGINI,NOMID LESOBL(1)='RLUX' LESOBL(2)='RLUY' LESOBL(3)='RLUZ' LESOBL(4)='RLRX' LESOBL(5)='RLRY' LESOBL(6)='RLRZ' LESOBL(7)='VX ' LESOBL(8)='VY ' LESOBL(9)='VZ ' * Elements de JOINTs GENE ELSE IF (MFR.EQ.55) THEN NBRFAC=1 SEGINI,NOMID LESFAC(1)='EPAI' * Macro element (element CIFL) ELSE IF (MFR.EQ.61)THEN NBROBL=2 SEGINI,NOMID LESOBL(1)='SECT' LESOBL(2)='INRZ' ENDIF * dans RIGI1.ESO : ajout de composantes facultatives pour le rendement NCARA=NBROBL NCARF=NBRFAC NCART=NCARA+NCARF MOCARA=NOMID MOTYPE = NOTYPE IF (MOCARA.NE.0) THEN CALL KOMCHA(IPCHE3,IPMAIL,CONM,MOCARA,MOTYPE,1,INFOS,3,IVACAR) SEGSUP,NOMID IF (IERR.NE.0) THEN KERRE=-ABS(IERR) GOTO 515 ENDIF IF (ISUPMA.EQ.1) THEN CALL VALCHE(IVACAR,NCART,IPMINT,IPPORE,MOCARA,MELE) IF (IERR.NE.0) THEN ISUPMA=0 KERRE=-ABS(IERR) GOTO 515 ENDIF ENDIF ENDIF IF (MOTYPE.NE.MOTYR8) SEGSUP,NOTYPE * dans RIGI1.ESO : 1) utilisation de la densite pour ponderer la prop * de phase si besoin, 2) MFR = 63, elements XFEM traites par RIGIXR *----------------------------------------------------------------------- *- 2.4 - Calcul de la matrice tangente selon le type d'element *----------------------------------------------------------------------- * 20 elements par ligne du GOTO GOTO (99,99,99, 4,99, 4,99, 4,99, 4,99,12,99, 4, 4, 4, 4,12,12,99, 1 99,22, 4, 4, 4, 4,27,28,29,30,99,99,99,99,35,35,35,35,35,35, 2 27,42,43,27,42,46,12,35,27,30,99,99,35,35,12,27,99,99,99,99, 3 99,99,99,99,99,99,99,99, 4, 4, 4, 4, 4,99,99,99,99,99,35,35, 4 35,35,35,84,85,86,42,42,99,99,99,42,27,12,46,42,42,42,99,99, 5 99,99,99,99,99,99,99,35,35,35,35,35,35,35,35,35,35,35,35,35, 6 35,35,42,42,42,42,42,99,99,99,99,99,99,99,99,99,99,99,99,99, 7 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,42,99,99,99, 8 99,99,99,99,99,99,99,42,42,42,42,42, 4, 4, 4, 4, 4, 4, 4, 4, 9 4, 4, 4, 4,99,99,99,99,99,99,99,99, 4, 4,99,99,99,99,99,99) & ,MELE IF (MELE.EQ.258.OR.MELE.EQ.260) GOTO 42 * Erreur : Element fini non encore implemente 99 CONTINUE KERRE=86 MOTERR(1:4)=NOMTP(MELE) MOTERR(5:12)='KTANGA ' GOTO 510 *----------------------------------------------------------------------- *-> Elements MASSIFs *----------------------------------------------------------------------- 4 CONTINUE NBNO=NBNN NBBB=NBNN SEGINI,WRK1,WRK2 IF (BMATE) THEN SEGINI,WTRAV NLG=NUMGEO(MELE) ENDIF IF (MAPL.EQ.5 .OR. MAPL.EQ.54) CALL ZDANUL(TRAC,LTRAC) IF (MAPL.EQ.54) SEGINI,WRK5 * Preparation a la recuperation de l'epaisseur MVALEP=0 DIM3=1.D0 IF (IFOUR.EQ.-2) THEN IF (IVACAR.NE.0) THEN MPTVAL=IVACAR MVALEP=IVAL(1) IF (MVALEP.GT.0) THEN MELVAL=MVALEP NELEP=VELCHE(/2) NPGEP=VELCHE(/1) ENDIF ENDIF ENDIF * Boucle sur les elements de la sous-zone ISOU DO 3004 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Calcul de la matrice de changement de repere (materiau a "tropie") IF (BMATE) THEN CALL RESHPT(1,NBNN,NLG,MELE,0,IPMIN2,IRT1) MINTE2=IPMIN2 SEGACT,MINTE2 NBSH=MINTE2.SHPTOT(/2) CALL RLOCAL(XE,MINTE2.SHPTOT,NBSH,NBNN,TXR) IF (NBSH.EQ.-1) THEN KERRE=525 GOTO 8004 ENDIF ENDIF * Mise a zero de la matrice de rigidite elementaire (IB) CALL ZERO(REL,LRE,LRE) * Cas des elements incompressibles : termes de la matrice B-BARRE IF (MFR.EQ.31) THEN CALL BBCALC(MELE,NBNN,MFR,IDIM,XE, & NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & NSTRS,LRE,IFOUR,NHRM,A,BB, & SHPTOT,SHPWRK,BGENE,XDPGE,YDPGE,PP,NOER) IF (NOER.NE.0) THEN CALL ERREUR(noer) RETURN ENDIF ENDIF * Boucle sur les points de Gauss de l'element IB ISDJC=0 DO 4004 IGAU=1,NBPGAU * Recuperation de l'epaisseur si donnee IF (MVALEP.GT.0) THEN MELVAL=MVALEP IBMN=MIN(IB ,NELEP) IGMN=MIN(IGAU,NPGEP) DIM3=VELCHE(IGMN,IBMN) ENDIF * Calcul de la matrice B et du jacobien DJAC au point de Gauss IGAU CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,DIM3, & XE,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) IF (DJAC.EQ.0.) THEN KERRE=259 INTERR(1)=IB GOTO 8004 ENDIF IF (DJAC.LT.0.) ISDJC=ISDJC+1 DJAC=ABS(DJAC)*POIGAU(IGAU) C En cas d'elements incompressibles : BGENE selon la methode B-BARRE IF (MFR.EQ.31) THEN CALL BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,NBNN,LRE,IFOUR,NSTRS,XE,DJAC,A,BB,BGENE) ENDIF IRET=0 * Recuperation des proprietes materielles utiles selon le modele MPTVAL=IVAMAT IF (MAPL.EQ.5) THEN MELVAL=IVAL(1) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) YYYY=VELCHE(IGMN,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB ,IELCHE(/2)) IGMN=MIN(IGAU,IELCHE(/1)) IMMM=IELCHE(IGMN,IBMN) CALL COTRA1(IMMM,YYYY,LTRAC,TRAC,NTRAC,KERRE) ELSE IF (MAPL.EQ.54) THEN MELVAL=IVAL(12) IBMN=MIN(IB ,IELCHE(/2)) IGMN=MIN(IGAU,IELCHE(/1)) ITREV=IELCHE(IGMN,IBMN) MELVAL=IVAL(13) IBMN=MIN(IB ,IELCHE(/2)) IGMN=MIN(IGAU,IELCHE(/1)) ICOEV=IELCHE(IGMN,IBMN) IPOS1=1 CALL COTRAB(ITREV,TRAC,LTRAC,IPOS1,0,NPOINT,KERRE) IF (KERRE.NE.0) THEN INTERR(1)=IB INTERR(2)=IGAU MOTERR(1:4)=NOMTP(MELE) GOTO 8004 ENDIF NTRAC1=NPOINT/2 IPOS2=IPOS1+NPOINT CALL COTRAB(ICOEV,TRAC,LTRAC,IPOS2,0,NPOINT,KERRE) NTRAC2=NPOINT/2 IRET=WRK5 ENDIF IF (KERRE.NE.0) THEN INTERR(1)=IB INTERR(2)=IGAU MOTERR(1:4)=NOMTP(MELE) GOTO 8004 ENDIF IF (BMATE) IRET=WTRAV * Contribution du pt de Gauss IGAU a la matrice tangente elementaire CALL DOHOT1(IVAMAT,NMATT,IVACON,NSTRS,IVARI,NVART,TRAC, & LTRAC,IGAU,IB,MATE,MAPL,XPREC,DTPS,IFOUR, & LHOOK,DDHOOK,IRET) IF (IRET.NE.1) THEN IF (IRET.EQ.-1) THEN KERRE=275 INTERR(1)=IB INTERR(2)=IGAU MOTERR(1:4)=NOMTP(MELE) C* ELSE IF (IRET.EQ.0) THEN ELSE KERRE=328 INTERR(1)=IFOUR MOTERR(1:8)=NOMAT(MATE) MOTERR(9:12)=NOMAC(MAPL) MOTERR(13:20)=NOMFR(MFR) ENDIF GOTO 8004 ENDIF IF (IRIGE7.EQ.2) THEN CALL BDBSTS(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ELSE CALL BDBST(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ENDIF 4004 CONTINUE * Fin de la boucle sur les points de Gauss IF (ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) THEN KERRE=195 INTERR(1)=IB GOTO 8004 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) IF (IRIGE7.EQ.2)THEN CALL REMPMS(REL,LRE,RE(1,1,IB)) ELSE CALL REMPMT(REL,LRE,RE(1,1,IB)) ENDIF 3004 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8004 CONTINUE * Menage local SEGSUP,WRK1,WRK2 IF (MAPL.EQ.54) SEGSUP,WRK5 IF (BMATE) SEGSUP,WTRAV GOTO 510 *----------------------------------------------------------------------- *-> Elements de raccord liquide-solide : * RAC2 LIA3 LIA4 RACO LICO LIC4 (MELE = 12 18 19 47 55 94) * => Elements SANS RIGIDITE (elastique & tangente) *----------------------------------------------------------------------- 12 CONTINUE C Les matrices elementaires sont nulles et ont ete mises a zero lors de C l'initialisation du segment XMATRI ! GOTO 510 *----------------------------------------------------------------------- *-> Element de type "Relations de conformites" (MELE=22) * Matrice TANGENTE = Matrice de RIGIDITE (ELASTIQUE) (cf. RIGI1) *----------------------------------------------------------------------- 22 CONTINUE IF (ITYPEL.NE.22) THEN KERRE=977 GOTO 510 ENDIF CALL RIGSUR(ISOU ,IPMATR, IMODEL) GOTO 510 *----------------------------------------------------------------------- *-> Element de type "Relations de conformites " (MELE=259) * Matrice TANGENTE = Matrice de RIGIDITE (ELASTIQUE) (cf. RIGI1) *----------------------------------------------------------------------- 259 CONTINUE IF (ITYPEL.NE.259) THEN KERRE=977 GOTO 510 ENDIF CALL RIGSUR(ISOU ,IPMATR, IMODEL) * Cas particulier si formulation X-FEM IF (IMODEL.INFELE(13).EQ.63) then CALL RIGSUX(ISOU ,IPMATR, IMODEL) ENDIF GOTO 510 *----------------------------------------------------------------------- *-> Elements COQ3 COQ8 COQ2 COQ4 COQ6 DST (MELE = 27 41 44 49 56 93) *-> Cas particulier : DKT elastique (MELE = 28 avec MAPL = 0) * Matrice TANGENTE = RIGIDITE ELASTIQUE (Appel a RIGI3) *----------------------------------------------------------------------- 27 CONTINUE MPTVAL=IVAMAT NBGMAT = 0 NELMAT = 0 IF (CMATE.EQ.'SECTION') THEN DO i=1,IVAL(/1) MELVAL=IVAL(i) IF (MELVAL.NE.0)THEN NBGMAT=MAX(NBGMAT,IELCHE(/1)) NELMAT=MAX(NELMAT,IELCHE(/2)) ENDIF ENDDO ELSE DO i=1,IVAL(/1) MELVAL=IVAL(i) IF (MELVAL.NE.0)THEN NBGMAT=MAX(NBGMAT,VELCHE(/1)) NELMAT=MAX(NELMAT,VELCHE(/2)) ENDIF ENDDO ENDIF CALL RIGI3(MATE,MELE,IPMAIL,IPMINT,IPMIN1,NBPGAU,LRE,NSTRS, & IVAMAT,IVACAR,CMATE,MFR,NBGMAT,NELMAT,1,LHOOK, & NMATT,LW,NPINT,IPMATR,IIPDPG) IF (IERR.NE.0) KERRE=-ABS(IERR) GOTO 510 *----------------------------------------------------------------------- *-> Element DKT (MELE = 28) *----------------------------------------------------------------------- 28 CONTINUE IF (MAPL.EQ.0) GOTO 27 NBNO=NBNN NBBB=NBNN SEGINI,WRK1,WRK2,WRK4 IF (MAPL.EQ.5) CALL ZDANUL(TRAC,LTRAC) * Preparation a la recuperation de caracteristiques geometriques MPTVAL=IVACAR EXCEN=0.D0 MVALEX=IVAL(2) IF (MVALEX.GT.0) THEN MELVAL=MVALEX NELEX=VELCHE(/2) NPGEX=VELCHE(/1) ENDIF CALFA=0.666666666666666666666666666666666666666666666D0 MVALCA=IVAL(NCART) IF (MVALCA.GT.0) THEN MELVAL=MVALCA NELCA=VELCHE(/2) NPGCA=VELCHE(/1) ENDIF * Boucle sur les elements de la sous-zone ISOU DO 3028 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Passage dans le repere local de l'element CALL VPAST(XE,BPSS) CALL VCORLC(XE,XEL,BPSS) * Mise a zero de la matrice de rigidite elementaire (IB) CALL ZERO(REL,LRE,LRE) * Boucle sur les points de Gauss de l'element IB DO 4028 IGAU=1,NBPGAU * Recuperation des caracteristiques geometriques MPTVAL=IVACAR MELVAL=IVAL(1) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) EPAIST=VELCHE(IGMN,IBMN) IF (MVALEX.GT.0) THEN MELVAL=MVALEX IBMN=MIN(IB ,NELEX) IGMN=MIN(IGAU,NPGEX) EXCEN=VELCHE(IGMN,IBMN) ENDIF IF (MVALCA.GT.0) THEN MELVAL=MVALCA IBMN=MIN(IB ,NELCA) IGMN=MIN(IGAU,NPGCA) CALFA=VELCHE(IGMN,IBMN) ENDIF * Calcul de la matrice B et du jacobien DJAC au point de Gauss IGAU CALL BMAT28(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,MFR,NBNO,LRE,IFOUR,NSTRS,0,1.D0,XEL, & SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE) DJAC=DJAC*POIGAU(IGAU) * Modification de la matrice B en cas d'excentrement non nul IF (ABS(EXCEN).GT.0.) THEN DO i=1,3 k=i+3 DO j=1,LRE BGENE(i,j)=BGENE(i,j)+EXCEN*BGENE(k,j) ENDDO ENDDO ENDIF * Recuperation des proprietes materielles utiles selon le modele IF (MAPL.EQ.5) THEN MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) YYYY=VELCHE(IGMN,IBMN) MELVAL=IVAL(3) IBMN=MIN(IB ,IELCHE(/2)) IGMN=MIN(IGAU,IELCHE(/1)) IMMM=IELCHE(IGMN,IBMN) CALL COTRA1(IMMM,YYYY,LTRAC,TRAC,NTRAC,KERRE) IF (KERRE.NE.0) THEN INTERR(1)=IB INTERR(2)=IGAU MOTERR(1:4)=NOMTP(MELE) GOTO 8028 ENDIF ENDIF * DOHOT3 se chargera de convertir les efforts generalises (IVACON) * et les variables internes generalisees (IVARI) en contraintes et * variables internes "locales" * Contribution du pt de Gauss IGAU a la matrice tangente elementaire CALL DOHOT3(IVAMAT,NMATT,IVACON,NSTRS,IVARI,NVART, & TRAC,LTRAC,CALFA,EPAIST,IGAU,IB,MATE,MAPL, & XPREC,DTPS,IFOUR,LHOOK,DDHOOK,IRET) IF (IRET.NE.1) THEN IF (IRET.EQ.-1) THEN KERRE=275 INTERR(1)=IB INTERR(2)=IGAU MOTERR(1:4)=NOMTP(MELE) C* ELSE IF (IRET.EQ.0) THEN ELSE KERRE=328 INTERR(1)=IFOUR MOTERR(1:8)=NOMAT(MATE) MOTERR(9:12)=NOMAC(MAPL) MOTERR(13:20)=NOMFR(MFR) ENDIF GOTO 8028 ENDIF CALL BDBST(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) 4028 CONTINUE * Fin de la boucle sur les points de Gauss * Calcul de la matrice tangente elementaire (IB) REL( 6, 6)=REL(5,5)*1.D-7 REL(12,12)=REL(6,6) REL(18,18)=REL(6,6) CALL TRANSK(REL,BPSS,LRE,3,1) * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3028 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8028 CONTINUE * Menage local SEGSUP,WRK1,WRK2,WRK4 GOTO 510 *----------------------------------------------------------------------- *-> Element POUTre (MELE=29) *-> Cas particulier : Element TIMO avec materiau ISOTROPE (MELE=84) *----------------------------------------------------------------------- 29 CONTINUE NBBB=NBNN SEGINI,WRK1,WRK3 * Boucle sur les elements de la sous-zone ISOU DO 3029 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Recuperation des caracteristiques geometriques (stockees dans WORK) MPTVAL=IVACAR DO IC=1,NCART r_z=0. IF (IVAL(IC).NE.0) THEN MELVAL=IVAL(IC) IBMN=MIN(IB,VELCHE(/2)) DO IGAU=1,NBNN IGMN=MIN(IGAU,VELCHE(/1)) r_z=r_z+VELCHE(IGMN,IBMN) ENDDO r_z=r_z/NBNN ENDIF WORK(IC)=r_z ENDDO * Recuperation des caracteristiques elastiques (stockees dans WORK) MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB ,VELCHE(/2)) YOUNG=VELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB ,VELCHE(/2)) XNU=VELCHE(1,IBMN) CISAIL=0.5*YOUNG/(1.+XNU) IF (BPLAN) THEN WORK(4)=YOUNG WORK(5)=CISAIL ELSE WORK(10)=YOUNG WORK(11)=CISAIL ENDIF * Modification de caracteristiques (INRY ou SECZ) selon les modeles MPTVAL=IVARI IF (MAPL.EQ.57.OR.MAPL.EQ.59) THEN MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) WORK(2)=VELCHE(1,IBMN)/YOUNG ELSE IF(MAPL.EQ.58.OR.MAPL.EQ.60) THEN MELVAL=IVAL(2) IBMN=MIN(IB,VELCHE(/2)) WORK(6)=VELCHE(1,IBMN)/CISAIL ENDIF * Calcul de la rigidite elementaire tangente (IB) IF (MELE.EQ.84) THEN IF (BPLAN) THEN CALL TIMRI2(REL,LRE,WORK,XE,WORK(12),IRET) ELSE CALL TIMRIG(REL,LRE,WORK,XE,WORK(12),IRET) ENDIF ELSE IF (BPLAN) THEN CALL POURI2(REL,LRE,WORK,XE,WORK(12),IRET) ELSE CALL POURIG(REL,LRE,WORK,XE,WORK(12),IRET) ENDIF ENDIF IF (IRET.NE.0) THEN IF (IRET.EQ.1) KERRE=128 IF (IRET.EQ.2) KERRE=138 INTERR(1)=ISOUS INTERR(2)=IB GOTO 8029 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3029 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8029 CONTINUE * Menage local SEGSUP,WRK1,WRK3 GOTO 510 *----------------------------------------------------------------------- *-> Elements linespring LISP et LISM en nonlineaire (MELE = 30 50) *----------------------------------------------------------------------- 30 CONTINUE IF (MAPL.EQ.0) GOTO 42 NBBB=NBNN SEGINI,WRK1,WRK3 * Boucle sur les elements de la sous-zone ISOU DO 3030 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) IE=1 * Recuperation des proprietes materielles (stockees dans WORK) IE1=IE MPTVAL=IVAMAT DO IC=1,NBPGAU DO i=1,2 MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO DO i=3,4 WORK(IE)=0.D0 IE=IE+1 ENDDO MELVAL=IVAL(3) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO * Recuperation des contraintes (stockees dans WORK) IE2=IE MPTVAL=IVACON DO IC=1,NBGS DO i=1,NSTRS MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB ,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO * Recuperation des caracteristiques geometriques (stockees dans WORK) IE3=IE MPTVAL=IVACAR DO IC=1,NBPGAU DO i=1,NCART MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO * Recuperation des variables internes (stockees dans WORK) IE4=IE MPTVAL=IVARI DO IC=1,NBGS DO i=1,NVART MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO IE5=IE+1 * Calcul de la rigidite elementaire tangente (IB) CALL LISPKT(XE,WORK(IE1),WORK(IE2),WORK(IE3),WORK(IE4),NSTRS, & NBGS,NBPGAU,MELE,XPREC,WORK(IE5),REL,I70,I157, & I158,IRET,KERRE) IF (IRET.EQ.-1.OR.KERRE.NE.0) THEN KERRE=270 IF (IRET.EQ.-1) KERRE=275 INTERR(1)=IB INTERR(2)=1 MOTERR(1:4)=NOMTP(MELE) GOTO 8030 ENDIF IF (I158.EQ.1) THEN KERRE=158 INTERR(1)=IB GOTO 8030 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3030 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8030 CONTINUE * Menage local SEGSUP,WRK1,WRK3 GOTO 510 *----------------------------------------------------------------------- *-> Elements liquide LTR3 LQU4 LCU8 LPR6 LTE4 LPY5 (MELE = 35 a 40) * Elements de surface libre LSU2 LSU3 LSU4 (MELE = 48 53 54) * Elements poreux TRIP QUAP CUBP TETP PRIP (MELE = 79 a 83) * Elements joints poreux JOP3 JOP6 JOP8 (MELE = 108 a 110) * Elements massifs polygonaux POLY (MELE = 111 a 122) * Matrice TANGENTE = RIGIDITE ELASTIQUE (Appel a RIGI2) *----------------------------------------------------------------------- 35 CONTINUE NCAR1 = NCART + 1 MPTVAL=IVAMAT NBGMAT = 0 NELMAT = 0 IF (CMATE.EQ.'SECTION') THEN DO i=1,IVAL(/1) IF (IVAL(i).NE.0) THEN MELVAL=IVAL(i) NBGMAT=MAX(NBGMAT,IELCHE(/1)) NELMAT=MAX(NELMAT,IELCHE(/2)) ENDIF ENDDO ELSE DO i=1,IVAL(/1) IF (IVAL(i).NE.0) THEN MELVAL=IVAL(i) NBGMAT=MAX(NBGMAT,VELCHE(/1)) NELMAT=MAX(NELMAT,VELCHE(/2)) ENDIF ENDDO ENDIF noer=0 CALL RIGI2 (MATE,MELE,IPMAIL,IPMINT,NBPGAU,LRE,NSTRS,IVAMAT, & IVACAR,CMATE,MFR,NBGMAT,NELMAT,1,LHOOK,NMATT,IPORE, & NDDL,IPMATR,IIPDPG,NCAR1,melpha,noer) IF (IERR.NE.0) KERRE=-ABS(IERR) GOTO 510 *----------------------------------------------------------------------- *-> Elements TUYA POI1 JOT3 JOI4 TRIH TUYO (MELE = 42 45 87 88 92 96) * LSE2 LITU BAR3 BAEX LIA2 QUAH CUBH (MELE = 97 98 123 a 127) * TRH6 JCT3 JCI4 JGI2 JGT3 JGI4 (MELE = 157 168 a 172) * CIFL SURE?? SHB8 (MELE = 258 259? 260) *-> Cas particuliers : LISP & LISM en elasticite (MELE= 30 50 & MAPL=0) * TUFI en elasticite (MELE=43 & MAPL=0) * Matrice TANGENTE = RIGIDITE ELASTIQUE (Appel a RIGI4) *----------------------------------------------------------------------- 42 CONTINUE MPTVAL=IVAMAT NBGMAT = 0 NELMAT = 0 IF (CMATE.EQ.'SECTION') THEN DO i=1,IVAL(/1) IF (IVAL(i).NE.0)THEN MELVAL=IVAL(i) NBGMAT=MAX(NBGMAT,IELCHE(/1)) NELMAT=MAX(NELMAT,IELCHE(/2)) ENDIF ENDDO ELSE DO i=1,IVAL(/1) IF (IVAL(i).NE.0)THEN MELVAL=IVAL(i) NBGMAT=MAX(NBGMAT,VELCHE(/1)) NELMAT=MAX(NELMAT,VELCHE(/2)) ENDIF ENDDO ENDIF CALL RIGI4(MATE,MELE,IPMAIL,IPMINT,NBPGAU,LRE,NSTRS, & IVAMAT,IVACAR,IVECT,CMATE,MFR,NBGMAT,NELMAT,1, & LHOOK,NMATT,NCART,ISOUS,LW,IPORE,IPMATR,IIPDPG) IF (IERR.NE.0) KERRE=-ABS(IERR) GOTO 510 *----------------------------------------------------------------------- *-> Element TUFI en non lineaire (MELE=43) *----------------------------------------------------------------------- 43 CONTINUE IF (MAPL.EQ.0) GOTO 42 C* Temporaire : si SIGF n'est pas definie, on utilise la matrice elastique MPTVAL=IVAMAT IF (IVAL(3).EQ.0) GOTO 42 NBBB=NBNN SEGINI,WRK1,WRK3 * Boucle sur les elements de la sous-zone ISOU DO 3043 IB=1,NBELEM * Recuperation des proprietes materielles (stockees dans WORK) IE =1 IE1=IE MPTVAL=IVAMAT DO IC=1,NBPGAU DO i=1,2 MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO DO i=3,4 WORK(IE)=0.D0 IE=IE+1 ENDDO MELVAL=IVAL(3) IF (MELVAL.NE.0) THEN IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) ELSE WORK(IE)=0.D0 ENDIF IE=IE+1 DO i=6,7 WORK(IE)=0.D0 IE=IE+1 ENDDO ENDDO * Recuperation des contraintes (stockees dans WORK) IE2=IE MPTVAL=IVACON DO IC=1,NBGS DO i=1,NSTRS MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO * Recuperation des caracteristiques geometriques (stockees dans WORK) IE3=IE MPTVAL=IVACAR DO IC=1,NBPGAU DO i=1,8 MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO WORK(IE)=0.D0 IE=IE+1 ENDDO * Recuperation des variables internes (stockees dans WORK) IE4=IE MPTVAL=IVARI DO IC=1,NBGS DO i=1,NVART MELVAL=IVAL(i) IGMN=MIN(IC,VELCHE(/1)) IBMN=MIN(IB,VELCHE(/2)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO ENDDO IE5=IE+1 * Calcul de la rigidite elementaire tangente (IB) CALL TUFIKT(WORK(IE1),WORK(IE2),WORK(IE3),WORK(IE4),REL,XPREC, & IRET) IF (IRET.NE.0) THEN INTERR(1)=ISOUS INTERR(2)=IB IF (IRET.EQ.1) KERRE=137 IF (IRET.EQ.2) KERRE=123 IF (IRET.EQ.3) KERRE=266 IF (IRET.EQ.4) THEN KERRE=275 INTERR(1)=IB INTERR(2)=1 MOTERR(1:4)=NOMTP(MELE) ENDIF GOTO 8043 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3043 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8043 CONTINUE * Menage local SEGSUP,WRK1,WRK3 GOTO 510 *----------------------------------------------------------------------- *-> Elements BARRe et CERCe (MELE = 46 95) *----------------------------------------------------------------------- 46 CONTINUE IF (MELE.EQ.95.AND.(IFOUR.NE.0.AND.IFOUR.NE.1)) GOTO 99 NBBB=NBNN SEGINI,WRK1,WRK3 MPTVAL=IVACAR MELVAL=IVAL(1) NELCAR=VELCHE(/2) MELCAR=IVAL(1) * Boucle sur les elements de la sous-zone ISOU DO 3046 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Recuperation de la section de l'element IB MELVAL=MELCAR IBMN=MIN(IB,NELCAR) SECT=VELCHE(1,IBMN) * Recuperation du module tangent selon le modele utilise MPTVAL=IVARI IF (MAPL.EQ.93) THEN MELVAL=IVAL(16) IBMN=MIN(IB,VELCHE(/2)) YOUNGT=VELCHE(1,IBMN) ELSE IF (MAPL.EQ.92) THEN MELVAL=IVAL(6) IBMN=MIN(IB,VELCHE(/2)) YOUNGT=VELCHE(1,IBMN) ELSE IF (MAPL.EQ.39) THEN MELVAL=IVAL(6) IBMN=MIN(IB,VELCHE(/2)) YOUNGT=VELCHE(1,IBMN) ELSE IF (MAPL.EQ.40) THEN MELVAL=IVAL(4) IBMN=MIN(IB,VELCHE(/2)) YOUNGT=VELCHE(1,IBMN) ELSE IF (MAPL.EQ.0) THEN MPTVAL=IVAMAT MELVAL=IVAL(1) IBMN=MIN(IB,VELCHE(/2)) YOUNGT=VELCHE(1,IBMN) ELSE KERRE=81 MOTERR(1:8)=CMATE MOTERR(9:16)=NOMFR(MFR/2+1) INTERR(1)=IFOUR GOTO 8046 ENDIF * Calcul de la rigidite elementaire tangente (IB) XHOOK=YOUNGT*SECT IF (MELE.EQ.46) THEN CALL BARRIG(REL,LRE,XHOOK,XE,IRET) IF (IRET.EQ.1) KERRE=128 ELSE IF (MELE.EQ.95) THEN CALL CERRIG(REL,LRE,XHOOK,XE,IRET) IF (IRET.EQ.1) KERRE=601 ENDIF IF (KERRE.NE.0) THEN INTERR(1)=ISOUS INTERR(2)=IB GOTO 8046 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3046 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8046 CONTINUE * Menage local SEGSUP,WRK1,WRK3 GOTO 510 *----------------------------------------------------------------------- *-> Element de poutre de TIMOschenko (MELE=84 - Formulation 'SECTION') *----------------------------------------------------------------------- 84 CONTINUE IF (CMATE.EQ.'ISOTROPE') GOTO 29 * Remarque : La formulation SECTION est le seul cas prevu actuellement. IF (CMATE.NE.'SECTION') THEN KERRE=193 MOTERR(1:8)=NOMFR(MFR) GOTO 510 ENDIF NBBB=NBNN SEGINI WRK1,WRK3 * Boucle sur les elements de la sous-zone ISOU DO 3084 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Recuperation des caracteristiques geometriques (stockees dans WORK) * Attention : on a tenu compte du fait que NCART=1 MPTVAL=IVACAR IF (IVAL(NCART).NE.0) THEN MELVAL=IVAL(NCART) IBMN=MIN(IB,IELCHE(/2)) IP=IELCHE(1,IBMN) IREF=(IP-1)*(IDIM+1) DO IC=1,IDIM WORK(IC)=XCOOR(IREF+IC) ENDDO ELSE DO IC=1,IDIM WORK(IC)=0.D0 ENDDO ENDIF MPTVAL=IVAMAT * Traitement dans le cas de la formulation section C** IF (CMATE.EQ.'SECTION') THEN * Dans le cas d'un materiau elastique, on prend les matrices de Hooke * si elles existent dans le MCHAML des proprietes materielles IF (MAPL.EQ.0.AND.IVAL(3).NE.0) THEN MELVAL=IVAL(3) IF (IB.LE.IELCHE(/2).OR.IELCHE(/1).GT.1) THEN IBMN=MIN(IB,IELCHE(/2)) MLREEL=IELCHE(1,IBMN) SEGACT,MLREEL CALL DOHOOO(PROG,LHOOK,DDHOOK) SEGDES,MLREEL ENDIF ELSE MELVAL=IVAL(1) IBMN=MIN(IB,IELCHE(/2)) IPMOD=IELCHE(1,IBMN) MELVAL=IVAL(2) IBMN=MIN(IB,IELCHE(/2)) IPCAR=IELCHE(1,IBMN) * Sinon calcul des matrices de Hooke a partir des proprietes elastiques IF (MAPL.EQ.0) THEN CALL FRIGIE(IPMOD,IPCAR,CRIGI,CMASS) * Ou calcul des matrices de Hooke a partir des variables internes ELSE MPTVAL=IVARI MELVAL=IVAL(2) IBMN=MIN(IB,IELCHE(/2)) IPVAR=IELCHE(1,IBMN) IF (IPVAR.NE.0) THEN CALL FRIGTA(IPMOD,IPCAR,IPVAR,CRIGI) ELSE CALL FRIGIE(IPMOD,IPCAR,CRIGI,CMASS) ENDIF ENDIF CALL DOHTIF(CRIGI,CMATE,IFOUR,LHOOK,DDHOOK,IRET) ENDIF IF (BPLAN) THEN CALL TIFRI2(REL,LRE,XE,WORK(12),LHOOK,DDHOOK,IRET) ELSE CALL TIFRIG(REL,LRE,WORK,XE,WORK(12),LHOOK,DDHOOK,IRET) ENDIF C** ENDIF IF (IRET.NE.0) THEN INTERR(1)=ISOUS INTERR(2)=IB IF (IRET.EQ.1) KERRE=128 IF (IRET.EQ.2) KERRE=138 GOTO 8084 ENDIF * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) CALL REMPMT(REL,LRE,RE(1,1,IB)) 3084 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8084 CONTINUE * Menage local SEGSUP,WRK1,WRK3 GOTO 510 *----------------------------------------------------------------------- *-> Element JOI2 - Materiau ISOTROPE (MELE = 85) *----------------------------------------------------------------------- *OF A voir : Erreur pour joi_ama.dgibi car MAPL=47 et on passe dans *OF DOUOTA qui ne traite que MAPL=35,56,91 85 CONTINUE IF (CMATE.NE.'ISOTROPE') THEN KERRE=834 MOTERR(1:8)=CMATE GOTO 510 ENDIF NBNO=NBNN NBBB=NBNN LW=100 SEGINI,WRK1,WRK2,WRK3,WRK4 * Boucle sur les elements de la sous-zone ISOU DO 3085 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Calcul des axes locaux de l'element IB CALL JO2LOC(XE,SHPTOT,NBNO,XEL,BPSS,NOQUAL) * Mise a zero de la matrice de rigidite elementaire (IB) CALL ZERO(REL,LRE,LRE) * Boucle sur les points de Gauss de l'element IB DO 4085 IGAU=1,NBPGAU * Calcul de la matrice B et du jacobien DJAC au point de Gauss IGAU CALL BJO2(IGAU,MFR,IFOUR,NIFOUR,XEL,BPSS,SHPTOT,SHPWRK, & BGENE,DJAC,IRET) DJAC=DJAC*POIGAU(IGAU) * Erreur si le jacobien est <= 0. IF (IRET.NE.0) THEN KERRE=612 INTERR(1)=IB GOTO 8085 ENDIF * Recuperation des proprietes materielles (stockees dans WORK) IE=1 IE1=IE MPTVAL=IVAMAT DO i=1,NMATT MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO * Calcul de la matrice de Hooke tangente au point de Gauss IGAU IF (MAPL.EQ.35.OR.MAPL.EQ.56.OR.MAPL.EQ.91) THEN * Recuperation des contraintes (stockees dans WORK) IE2=IE MPTVAL=IVACON DO i=1,NSTRS MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO * Recuperation des variables internes (stockees dans WORK) IE3=IE MPTVAL=IVARI DO i=1,NVARI+NVARF MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO CALL DOUOTA(WORK(IE1),CMATE,IFOUR,MAPL,WORK(IE2), & WORK(IE3),LHOOK,DDHOOK,IRET) ELSE CALL DOUO88(WORK(IE1),CMATE,IFOUR,LHOOK,DDHOOK,IRET) ENDIF IF (IRET.EQ.0) THEN KERRE=81 INTERR(1)=IFOUR MOTERR(1:8)=CMATE MOTERR(9:16)=NOMFR(MFR/2+1) GOTO 8085 ENDIF * Contribution du pt de Gauss IGAU a la matrice tangente elementaire IF (IRIGE7.EQ.2) THEN CALL BDBSTS(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ELSE CALL BDBST(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ENDIF 4085 CONTINUE * Fin de la boucle sur les points de Gauss * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) IF (IRIGE7.EQ.2) THEN CALL REMPMS(REL,LRE,RE(1,1,IB)) ELSE CALL REMPMT(REL,LRE,RE(1,1,IB)) ENDIF 3085 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8085 CONTINUE * Menage local SEGSUP,WRK1,WRK2,WRK3,WRK4 GOTO 510 *----------------------------------------------------------------------- *-> Element JOI3 - Materiau ISOTROPE (MELE = 86) *----------------------------------------------------------------------- 86 CONTINUE IF (CMATE.NE.'ISOTROPE') THEN KERRE=834 MOTERR(1:8)=CMATE GOTO 510 ENDIF NBNO=NBNN NBBB=NBNN LW=100 SEGINI,WRK1,WRK2,WRK3,WRK4 * Boucle sur les elements de la sous-zone ISOU DO 3086 IB=1,NBELEM * Recuperation des coordonnees des noeuds de l'element IB CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) * Mise a zero de la matrice de rigidite elementaire (IB) CALL ZERO(REL,LRE,LRE) * Boucle sur les points de Gauss de l'element IB DO 4086 IGAU=1,NBPGAU * Calcul de la matrice B et du jacobien DJAC au point de Gauss IGAU CALL JO3LOC(XE,SHPTOT,IGAU,NBNO,BPSS) CALL BJO3(IGAU,MFR,IFOUR,NIFOUR,XE,BPSS,SHPTOT,SHPWRK, & BGENE,DJAC,IRET) DJAC=DJAC*POIGAU(IGAU) * Erreur si le jacobien est <= 0. IF (IRET.NE.0) THEN KERRE=612 INTERR(1)=IB CALL ERREUR(612) GOTO 8086 ENDIF * Recuperation des proprietes materielles (stockees dans WORK) IE=1 IE1=IE MPTVAL=IVAMAT DO i=1,NMATT MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO * Calcul de la matrice de Hooke tangente au point de Gauss IGAU IF (MAPL.EQ.35.OR.MAPL.EQ.56.OR.MAPL.EQ.91) THEN * Recuperation des contraintes (stockees dans WORK) IE2=IE MPTVAL=IVACON DO i=1,NSTRS MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO * Recuperation des variables internes (stockees dans WORK) IE3=IE MPTVAL=IVARI DO i=1,NVARI+NVARF MELVAL=IVAL(i) IBMN=MIN(IB ,VELCHE(/2)) IGMN=MIN(IGAU,VELCHE(/1)) WORK(IE)=VELCHE(IGMN,IBMN) IE=IE+1 ENDDO CALL DOUOTA(WORK(IE1),CMATE,IFOUR,MAPL,WORK(IE2), & WORK(IE3),LHOOK,DDHOOK,IRET) ELSE CALL DOUO88(WORK(IE1),CMATE,IFOUR,LHOOK,DDHOOK,IRET) ENDIF IF (IRET.EQ.0) THEN KERRE=81 INTERR(1)=IFOUR MOTERR(1:8)=CMATE MOTERR(9:16)=NOMFR(MFR/2+1) GOTO 8086 ENDIF * Contribution du pt de Gauss IGAU a la matrice tangente elementaire IF (IRIGE7.EQ.2) THEN CALL BDBSTS(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ELSE CALL BDBST(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL) ENDIF 4086 CONTINUE * Remplissage de la matrice de rigidite elementaire (RE(.,.,IB)) IF (IRIGE7.EQ.2) THEN CALL REMPMS(REL,LRE,RE(1,1,IB)) ELSE CALL REMPMT(REL,LRE,RE(1,1,IB)) ENDIF 3086 CONTINUE * Fin de la boucle sur les elements * Etiquette de gestion des erreurs 8086 CONTINUE * Menage local SEGSUP,WRK1,WRK2,WRK3,WRK4 GOTO 510 *----------------------------------------------------------------------- * desactivation des segments propres a la zone geometrique ISOU *----------------------------------------------------------------------- 510 CONTINUE IF (IPMATR.NE.0) THEN * Symetrisation de la matrice si demandee IF (IRIGE7.EQ.2 .AND. IKTSYM.EQ.1) THEN LRE = RE(/1) DO 5100 IB=1,NBELEM DO i=1,LRE DO j=i+1,LRE RE(i,j,IB) = 0.5*( RE(i,j,IB)+RE(j,i,IB) ) RE(j,i,IB) = RE(i,j,IB) ENDDO ENDDO 5100 CONTINUE ENDIF ENDIF IF (IPDSCR.NE.0) SEGDES,DESCR 515 CONTINUE IF (ISUPMA.EQ.1) THEN CALL DTMVAL(IVACAR,3) ELSE CALL DTMVAL(IVACAR,1) ENDIF IF (ISUPMA.EQ.1) THEN CALL DTMVAL(IVAMAT,3) ELSE CALL DTMVAL(IVAMAT,1) ENDIF IF (ISUPVA.EQ.1) THEN CALL DTMVAL(IVARI,3) ELSE CALL DTMVAL(IVARI,1) ENDIF IF (ISUPCO.EQ.1) THEN CALL DTMVAL(IVACON,3) ELSE CALL DTMVAL(IVACON,1) ENDIF ** MELEME=IPMAIL * Mise a jour du segment MRIGID en cas de SUCCES ! IF (KERRE.EQ.0) THEN COERIG(ISOU)=1.D0 IRIGEL(1,ISOU)=IPMAGD IRIGEL(2,ISOU)=0 IRIGEL(3,ISOU)=IPDSCR IRIGEL(4,ISOU)=IPMATR IRIGEL(5,ISOU)=NIFOUR IRIGEL(6,ISOU)=0 IRIGEL(7,ISOU)=IRIGE7*(1-IKTSYM) xmatri.symre=irigel(7,isou) SEGDES,XMATRI IRIGEL(8,ISOU)=0 * Sinon en cas d'ERREUR : sortie prematuree ! ELSE * Affichage erreur si KERRE > 0 (car KERRE<0 qd erreur deja imprimee) IF (KERRE.GT.0) CALL ERREUR(KERRE) IF (IPDSCR.NE.0) SEGSUP,DESCR IF (IPMATR.NE.0) SEGSUP,XMATRI ** IF (IPMAGD.GT.0.AND.IPMAGD.NE.IPMAIL) THEN ** MELEME=IPMAGD ** SEGSUP,MELEME ** ENDIF GOTO 551 ENDIF * Fin de la boucle (5000) de PARTITIONNEMENT du segment XMATRI 5000 CONTINUE *----------------------------------------------------------------------- * Fin de la boucle sur les sous-zones du modele *----------------------------------------------------------------------- 500 CONTINUE * Sortie du sous-programme 551 CONTINUE IF (KERRE.EQ.0) THEN SEGDES,MRIGID IPRIGI=MRIGID ELSE SEGSUP,MRIGID IPRIGI=0 ENDIF * Desactivation du modele "deroule" 550 CONTINUE MMODEL = IPMODL SEGDES,MMODEL meleme = MAILDG IF (meleme.NE.0) SEGDES,meleme c RETURN END