C MENISM SOURCE CB215821 20/11/25 13:34:22 10792 SUBROUTINE MENISM ********************************************************************* * * T_CHPO1 T_CHML1=MESM MOD1 BLOQ1 (('TOUT') or 'ROTA'); * * CALCUL DES MECHANISMES ELEMENTAIRES D'UNE STRUCTURE * COMPOSEE DE POUTRE. * * P.PEGON (ISPRA) AOUT 1996 ********************************************************************* IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMRIGID -INC SMMODEL -INC SMCOORD -INC SMELEME -INC SMCHAML -INC SMCHPOI -INC SMTABLE SEGMENT ICPR(nbpts) SEGMENT WORK0 REAL*8 C(MRIGI,2*NCOMP) REAL*8 Q(2*NCOMP,2*NCOMP) REAL*8 Q1(2*NCOMP,2*NCOMP+MECAN) REAL*8 C1(MRIGI,2*NCOMP) REAL*8 C2(MRIGI,2*NCOMP+MECAN) REAL*8 XE(3,2) REAL*8 V1(2*NCOMP) REAL*8 TH(MRIGI) ENDSEGMENT SEGMENT WORK1 REAL*8 A(NEQUA,NINCO) ENDSEGMENT SEGMENT WORK2 REAL*8 B(NINCO,NINCO) ENDSEGMENT SEGMENT WORK3 REAL*8 B1(NNNCO,NMECA) ENDSEGMENT * PARAMETER(NCOMP=3,NECAN=3,MRIGI=3) CHARACTER*(LOCOMP) COMP(NCOMP),CCOMP,CANI(NECAN) PARAMETER (NMOCLE=2) CHARACTER*4 MOTCLE(NMOCLE) * DATA COMP/'UX ','UY ','RZ '/ DATA CANI/'RZP1','RZP2','UP2 '/ DATA MOTCLE/'TOUT','ROTA'/ * * Lecture des donnees * CALL LIROBJ('MMODEL ',MMODEL,1,IRET) IF(IRET.EQ.0)RETURN CALL LIROBJ('RIGIDITE',MRIGID,1,IRET) IF(IRET.EQ.0)RETURN CALL LIRMOT(MOTCLE,NMOCLE,IVAL,0) IF(IVAL.EQ.0.OR.IVAL.EQ.1) THEN MECAN=3 ELSE IF(IVAL.EQ.2)THEN MECAN=2 ENDIF * * On verifie que l'on est en 2D et que les elements sont de * poutre ou de timo * IF (IDIM.NE.2)THEN WRITE(IOIMP,*)'MENISM: ne fonctionne que en 2D' CALL ERREUR(19) RETURN ENDIF * IRET=0 SEGACT,MMODEL NZONE=KMODEL(/1) DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) SEGACT,IMODEL IF(NEFMOD.NE.29.AND.NEFMOD.NE.84)THEN IRET=IRET+1 ENDIF ENDDO IF(IRET.GT.0)THEN WRITE(IOIMP,*) > 'MENISM: ne fonctionne que pour des POUT ou des TIMO' CALL ERREUR(16) GOTO 9999 ENDIF * * On numerote les points support actifs * SEGINI,ICPR DO I=1,nbpts ICPR(I)=0 ENDDO NPOI=0 NBEL=0 DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) MELEME=IMAMOD SEGACT,MELEME NBELEM=NUM(/2) NBEL=NBEL+NBELEM DO I=1,NUM(/1) DO J=1,NBELEM IKI=NUM(I,J) IF (ICPR(IKI).EQ.0)THEN NPOI=NPOI+1 ICPR(IKI)=NPOI ENDIF ENDDO ENDDO ENDDO * * Verification que les c.l. s'appliquent a des points * et a des d.o.f.s du probleme * NLIM=0 SEGACT,MRIGID NRI=IRIGEL(/1) NR=IRIGEL(/2) DO I=1,NR MELEME=IRIGEL(1,I) DESCR =IRIGEL(3,I) xMATRI=IRIGEL(4,I) NEGALI=IRIGEL(6,I) IF (NEGALI.NE.0)THEN WRITE(IOIMP,*) > 'MENISM: on ne veut que des rigidites d"egalites' CALL ERREUR(21) GOTO 9997 ENDIF SEGACT,MELEME IF (ITYPEL.NE.22)THEN WRITE(IOIMP,*) > 'MENISM: on ne veut que des rigidites de blocage' CALL ERREUR(21) GOTO 9997 ENDIF NBELEM=NUM(/2) NBNN=NUM(/1) NLIM=NLIM+NBELEM DO J=1,NBELEM DO K=2,NBNN IF(ICPR(NUM(K,J)).EQ.0)THEN WRITE(IOIMP,*)'MENISM: ', > 'un point soumis a blocage n"est pas dans le modele' CALL ERREUR(21) GOTO 9997 ENDIF ENDDO ENDDO SEGACT,xMATRI NELRIG=re(/3) IF(NBELEM.NE.NELRIG)THEN WRITE(IOIMP,*) > 'MENISM: il y a un probleme sioux dans une des rigidites' CALL ERREUR(21) GOTO 9997 ENDIF SEGACT,DESCR NLIGRP=LISINC(/2) IF(NBNN.NE.NLIGRP)THEN WRITE(IOIMP,*) > 'MENISM: il y a un probleme sioux dans une des rigidites' CALL ERREUR(21) GOTO 9997 ENDIF DO K=2,NLIGRP CCOMP=LISINC(K) DO 1 J=1,NCOMP IF(COMP(J).EQ.CCOMP)GOTO 2 1 CONTINUE WRITE(IOIMP,*) > 'MENISM: un des d.o.f. de blocage n"est pas dans le modele' CALL ERREUR(21) GOTO 9997 2 CONTINUE ENDDO ENDDO * * On est en mesure de commencer le travail en se souvenant que * * NCOMP=nb de d.o.f par point * MECAN=nb de mecanisme par element * MRIGI=nb de contrainte rigide * NPOI=nb de point support de NCOMP dof * NBEL=nb d'elements * NLIM=nb de condition aux limites * * En consequence on peut allouer le pb avec * * NINCO=nb d'inconnues * NEQUA=nb d'equations * NINCO=NCOMP*NPOI+MECAN*NBEL NEQUA=MRIGI*NBEL+NLIM NNNCO=NCOMP*NPOI NMECA=NINCO-NEQUA * * On commence par la matrice A * SEGINI,WORK1 CALL ZERO(A,NEQUA,NINCO) * * Boucle sur les elements * SEGINI,WORK0 IEL=0 DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) MELEME=IMAMOD DO I=1,NUM(/2) IEL=IEL+1 * * Calcul de C2 elementaire * CALL DOXE(XCOOR,IDIM,2,NUM,I,XE) CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP) CALL TIPOQ1(Q,2*NCOMP, Q1,2*NCOMP+MECAN) CALL MATMAT(C,Q1,MRIGI,2*NCOMP,2*NCOMP+MECAN,C2) * * Assemblage de C2 dans A * * 1) noeud 1 de l'element * ILIGN=MRIGI*(IEL-1) ICOLO=NCOMP*(ICPR(NUM(1,I))-1) DO J=1,MRIGI DO K=1,NCOMP A(ILIGN+J,ICOLO+K)=C2(J,K) ENDDO ENDDO * * 2) noeud 2 de l'element * ICOLO=NCOMP*(ICPR(NUM(2,I))-1) DO J=1,MRIGI DO K=1,NCOMP A(ILIGN+J,ICOLO+K)=C2(J,K+NCOMP) ENDDO ENDDO * * 3) degres suplementaires de l'element * ICOLO=NCOMP*NPOI+MECAN*(IEL-1) DO J=1,MRIGI DO K=1,MECAN A(ILIGN+J,ICOLO+K)=C2(J,K+2*NCOMP) ENDDO ENDDO ENDDO ENDDO * * Assemblage des conditions aux limites * ILIGN=MRIGI*NBEL NRI=IRIGEL(/1) NR=IRIGEL(/2) DO I=1,NR MELEME=IRIGEL(1,I) DESCR=IRIGEL(3,I) xMATRI=IRIGEL(4,I) DO J=1,NUM(/2) ILIGN=ILIGN+1 * XMATRI=IMATTT(J) * SEGACT,XMATRI DO K=2,NUM(/1) CCOMP=LISINC(K) DO 3 L=1,NCOMP IF(COMP(L).EQ.CCOMP)GOTO 4 3 CONTINUE 4 CONTINUE ICOLO=NCOMP*(ICPR(NUM(K,J))-1) A(ILIGN,ICOLO+L)=RE(1,K,j) ENDDO * SEGDES,XMATRI ENDDO ENDDO * * On alloue B et on appelle le solveur * SEGINI,WORK2 CALL GAJOME(A,NEQUA,NINCO,B,IRET,IOIMP) SEGSUP,WORK1 IF (IRET.NE.0)THEN SEGSUP,WORK2 WRITE(IOIMP,*) > 'MENISM: il y a un mecanisme global' CALL ERREUR(26) GOTO 9997 ENDIF * * On ne garde que la partie active de B dans B1 * SEGINI,WORK3 DO J=1,NMECA DO I=1,NNNCO B1(I,J)=B(I,J+NEQUA) ENDDO ENDDO SEGSUP,WORK2 * * On prepare les tables de sortie * M=NMECA+1 SEGINI,MTABLE MLOTAB=M MTABTI(1)='MOT' CALL POSCHA('SOUSTYPE',IRET) MTABII(1)=IRET MTABTV(1)='MOT' CALL POSCHA('MECANISMES_PAR_NOEUDS',IRET) MTABIV(1)=IRET DO I=1,NMECA MTABTI(I+1)='ENTIER' MTABII(I+1)=I MTABTV(I+1)='CHPOINT ' MTABIV(I+1)=0 ENDDO MTAB1=MTABLE * SEGINI,MTABLE=MTAB1 CALL POSCHA('MECANISMES_PAR_ELEMENTS',IRET) MTABIV(1)=IRET DO I=1,NMECA MTABTV(I+1)='MCHAML ' ENDDO MTAB2=MTABLE * * Maillage support des chanpoint solution * NBSOUS=0 NBREF=0 NBELEM=NPOI NBNN=1 SEGINI IPT1 IPT1.ITYPEL=1 DO I=1,nbpts IF (ICPR(I).NE.0)THEN IPT1.NUM(1,ICPR(I))=I ENDIF ENDDO SEGDES,IPT1 * * Generation des chanpoint solution * DO I=1,NMECA IF(I.EQ.1)THEN * prototype MCHPOI NAT=2 NSOUPO=1 SEGINI,MCHPOI MOCHDE=' CHPOINT CREE PAR MENISM' MTYPOI='MECANISM' JATTRI(1)=1 IFOPOI=IFOUR MCHPO1=MCHPOI * prototype MSOUPO NC=NCOMP SEGINI,MSOUPO DO J=1,NCOMP NOCOMP(J)=COMP(J) ENDDO IGEOC=IPT1 MSOUP1=MSOUPO ELSE SEGINI,MCHPOI=MCHPO1 SEGINI,MSOUPO=MSOUP1 ENDIF IPCHP(1)=MSOUPO SEGDES,MCHPOI N=NPOI SEGINI,MPOVAL IPOVAL=MPOVAL SEGDES,MSOUPO DO J=1,NPOI DO K=1,NCOMP VPOCHA(J,K)=B1(NCOMP*(J-1)+K,I) ENDDO ENDDO SEGDES,MPOVAL MTAB1.MTABIV(I+1)=MCHPOI ENDDO SEGDES,MTAB1 * * Generation des mchaml solution * DO I=1,NMECA * * Chapeau * IF(I.EQ.1)THEN * prototype MCHAML N2=MECAN SEGINI,MCHAML DO J=1,MECAN NOMCHE(J)=CANI(J) TYPCHE(J)='REAL*8 ' ENDDO MCHAM1=MCHAML * prototype MCHELM N1=NZONE L1=9 N3=6 SEGINI,MCHELM TITCHE='MECANISME' IFOCHE=IFOUR DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) IMACHE(IZONE)=IMAMOD INFCHE(IZONE,1)=1 INFCHE(IZONE,2)=0 INFCHE(IZONE,3)=NIFOUR INFCHE(IZONE,4)=0 INFCHE(IZONE,5)=0 INFCHE(IZONE,6)=0 ENDDO MCHEL1=MCHELM * N1PTEL=1 N2PTEL=0 N2EL=0 ELSE SEGINI,MCHELM=MCHEL1 ENDIF * * Boucle sur les elements * DO IZONE=1,NZONE SEGINI,MCHAML=MCHAM1 ICHAML(IZONE)=MCHAML IMODEL=KMODEL(IZONE) MELEME=IMAMOD N1EL=NUM(/2) DO J=1,MECAN SEGINI,MELVAL IELVAL(J)=MELVAL ENDDO DO J=1,N1EL * * Calcul de C1 elementaire * CALL DOXE(XCOOR,IDIM,2,NUM,J,XE) CALL TIPOCQ(XE, C,Q,MRIGI,2*NCOMP) CALL MATMAT(C,Q,MRIGI,2*NCOMP,2*NCOMP,C1) * * On remplit V1 * * 1) noeud 1 de l'element * ICOLO=NCOMP*(ICPR(NUM(1,J))-1) DO K=1,NCOMP V1(K)=B1(ICOLO+K,I) ENDDO * * 2) noeud 2 de l'element * ICOLO=NCOMP*(ICPR(NUM(2,J))-1) DO K=1,NCOMP V1(K+NCOMP)=B1(ICOLO+K,I) ENDDO * * Calcul de THeta elementaire * CALL MATVE1(C1,V1,MRIGI,2*NCOMP,TH,2) * * Remplissage des MELVAL * DO K=1,MECAN MELVAL=IELVAL(K) VELCHE(1,J)=TH(K) ENDDO ENDDO DO J=1,MECAN MELVAL=IELVAL(J) SEGDES,MELVAL ENDDO SEGDES,MCHAML ENDDO SEGDES,MCHELM MTAB2.MTABIV(I+1)=MCHELM ENDDO SEGDES,MTAB2 SEGSUP,MCHAM1 * * A faire ... * * WARNING!!!!!!!!! SIGNE DES ROTATIONS!!!!!!! fait en modifiant C SEGSUP,WORK3 SEGSUP,WORK0 * * On rend la main a GIBIANE * CALL ECROBJ('TABLE',MTAB2) CALL ECROBJ('TABLE',MTAB1) * * Desactivation de la rigidite (imatri, descripteur et maillage) * 9997 DO I=1,NR MELEME=IRIGEL(1,I) SEGDES,MELEME DESCR=IRIGEL(3,I) SEGDES,DESCR xMATRI=IRIGEL(4,I) SEGDES,xMATRI ENDDO SEGDES,MRIGID * * Supression de ICPR * SEGSUP,ICPR * * Desactivation des maillages du modele * DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) MELEME=IMAMOD SEGDES,MELEME ENDDO * * Desactivation du modele * 9999 DO IZONE=1,NZONE IMODEL=KMODEL(IZONE) SEGDES,IMODEL ENDDO SEGDES,MMODEL * RETURN END