bloque
C BLOQUE SOURCE FANDEUR 22/01/03 21:15:03 11237 C----------------------------------------------------------------------- C Cet operateur impose les BLOCAGES C C Syntaxe : C C ENC1 = BLOQUER ( DEPL ) ( ROTA ) POI1 C C ou ENC1 = BLOQUER RADIAL P1 (P2) MELEME C ORTHOR P1 (P2) MELEME C C ou ENC1 = BLOQUER (DEPL) (ROTA) DIRECTION V1 MELEME C C DIM = 1 ( UX UY UZ ) ou ( UR UZ ) | C DIM = 2 OU 3 ( UX UY UZ RX RY RZ ) | MELEME C AXISYM ( RX RZ RT UT ) | Clist t C POI1 = OBJET DE TYPE POINT (entree) C MELEME = OBJET DE TYPE MELEME (entree) C ENC1 = OBJET DE TYPE RIGIDITE (sortie) C C 1) On peut imposer des BLOQUAGES UNILATERAUX en specifiant les C mots-cles MINIMUM ou MAXIMUM. C 2) On peut imposer des BLOQUAGES de type FROTTEMENT en specifiant C le mot-cle FROTTEMENT. C----------------------------------------------------------------------- C Juillet 2003 : passage a un seul multiplicateur C----------------------------------------------------------------------- SUBROUTINE BLOQUE IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCGEOME -INC CCREEL -INC CCHAMP -INC SMCHPOI -INC SMCOORD -INC SMELEME -INC SMLMOTS POINTEUR LMDEPL.MLMOTS -INC SMMODEL -INC SMRIGID -INC SMTABLE -INC TMTRAV SEGMENT MSWBLO CHARACTER*4 MOTDDL(0) ENDSEGMENT segment lispoi INTEGER pilpoi(mpoin),pilmul(mpoin) endsegment DIMENSION XNOR(3),U1(3),U2(3) CHARACTER*4 CHADDL CHARACTER*8 MOTRIG CHARACTER*4 MOTPV(3), MOTBLO(5) CHARACTER*4 MODEPL(6),MODEDU(6), MOROTA(5),MORODU(5) CHARACTER*4 MODE1D(2),MOFO1D(2) DATA MOTRIG / 'RIGIDITE' / DATA MOTPV / 'MINI','MAXI','FROT' / DATA MOTBLO / 'DEPL','ROTA','RADI','ORTH','DIRE' / DATA MODEPL / 'UX ','UY ','UZ ','UR ','UZ ','UT ' / DATA MODEDU / 'FX ','FY ','FZ ','FR ','FZ ','FT ' / DATA MOROTA / 'RX ','RY ','RZ ','RT ','RS ' / DATA MORODU / 'MX ','MY ','MZ ','MT ','MS ' / C Tableaux MODE1D et MOFO1D sont utilises pour certains modes 1D DATA MODE1D / 'UX ','UZ ' / DATA MOFO1D / 'FX ','FZ ' / C Pour ne pas avoir de verrouillage sur MCOORD en // SEGDES,MCOORD SEGACT,MCOORD*MOD c lecture table liaisons statiques c traitement table liaisons statiques IF (iretou.NE.0) THEN RETURN ENDIF C Est-ce une condition unilaterale ? NILATE=0 IF (IPO.EQ.1) NILATE=-1 IF (IPO.EQ.2) NILATE=1 IF (IPO.EQ.3) NILATE=2 C Pas de frottement en 1D IF (IPO.EQ.3.AND.IDIM.EQ.1) THEN INTERR(1)=IDIM MOTERR(1:4)=MOTPV(3) RETURN ENDIF idimp1=IDIM+1 C Quelques initialisations selon le type de probleme C Cas IDIM = 1 : C ISPE1D = 1 si IDIM=1 et IFOUR=9 ou 10, car les noms de DDL primaux et C variables duales ne sont pas dans l'ordre "classique" (un traitement C specifique est alors necessaire). ISPE1D=0 C Deformations planes ou contraintes planes ou defo. plane gene : IF (IFOUR.EQ.-1.OR.IFOUR.EQ.-2.OR.IFOUR.EQ.-3) THEN LDEPL=2 IADEPL=0 LROTA=1 IAROTA=2 C Axisymetrique : ELSE IF (IFOUR.EQ.0) THEN LDEPL=2 IADEPL=3 LROTA=1 IAROTA=3 C Fourier : ELSE IF (IFOUR.EQ.1) THEN LDEPL=3 IADEPL=3 LROTA=1 IAROTA=3 C Tridimensionnel : ELSE IF (IFOUR.EQ.2) THEN LDEPL=3 IADEPL=0 LROTA=3 IAROTA=0 C Massif 1D (IDIM=1) : ELSE IF (IFOUR.GE.3.AND.IFOUR.LE.15) THEN IF (IFOUR.LE.6) THEN LDEPL=1 IADEPL=0 ELSE IF (IFOUR.GE.7.AND.IFOUR.LE.10) THEN LDEPL=2 IADEPL=0 IF (IFOUR.EQ.9.OR.IFOUR.EQ.10) ISPE1D=1 ELSE IF (IFOUR.EQ.11) THEN LDEPL=3 IADEPL=0 ELSE IF (IFOUR.EQ.15) THEN LDEPL=2 IADEPL=3 ELSE LDEPL=1 IADEPL=3 ENDIF LROTA=0 IAROTA=0 C Autres cas : ELSE LDEPL=0 IADEPL=0 LROTA=0 IAROTA=0 ENDIF C Initialisation de la liste des ddls MOTDDL (segment MSWBLO) SEGINI,MSWBLO C----------------------------------------------------------------------- C Lecture eventuelle des MOTS autres que des DDLS C----------------------------------------------------------------------- C Lecture eventuelle de 'DEPL' et/ou 'ROTA' C -------------------- IDEPL=0 IROTA=0 IF (IMOT.EQ.1) IDEPL=1 IF (IMOT.EQ.2) IROTA=1 IF (IMOT.NE.0) GOTO 481 C Lecture eventuelle de 'RADI','ORTH' ou 'DIRE' C -------------------- IRADIA=0 IDIREC=0 IDICHP=0 IF (IMOT.EQ.0) THEN IF (IDEPL.EQ.1) GOTO 44801 IF (IROTA.EQ.1) GOTO 44802 IBDDL=MOTDDL(/2) IF (IBDDL.NE.0) GOTO 449 IF (IBDDL.EQ.0) GOTO 445 ENDIF C En DIMENSION 1, les mots-cles 'RADI,'ORTH' et 'DIRE' sont interdits. IF (IDIM.EQ.1) THEN INTERR(1)=IDIM MOTERR(1:4)=MOTBLO(2+IMOT) GOTO 1000 ENDIF GO TO (44803,44803,44804),IMOT C Traitement des mots-cles : Mise a jour de MOTDDL C ---------------------------- C On a trouve le mot DEPLAcement 44801 IDEPL=0 C Cas particulier pour certains modes de IDIM=1 IF (ISPE1D.EQ.1) THEN DO i=1,LDEPL MOTDDL(**)=MODE1D(IADEPL+i) MOTDDL(**)=MOFO1D(IADEPL+i) ENDDO C Cas general ELSE DO i=1,LDEPL MOTDDL(**)=MODEPL(IADEPL+i) MOTDDL(**)=MODEDU(IADEPL+i) ENDDO ENDIF GOTO 4480 C On a trouve le mot ROTAtion 44802 IROTA=0 DO i=1,LROTA MOTDDL(**)=MOROTA(IAROTA+i) MOTDDL(**)=MORODU(IAROTA+i) ENDDO GOTO 4480 C On a trouve le mot RADial ou le mot ORTHoradial 44803 IRADIA=IMOT IF (IRETOU.EQ.0) GOTO 1000 j=(KPOINT-1)*idimp1 DO i=1,IDIM U1(i)=XCOOR(j+i) ENDDO IF (IDIM.EQ.3) THEN C Lecture du 2nd point de l'axe (en 3D) IF (IRETOU.EQ.0) GOTO 1000 j=(KPOINT-1)*idimp1 YL=0.D0 DO i=1,IDIM U2(i)=XCOOR(j+i)-U1(i) YL=YL+U2(i)*U2(i) ENDDO C Calcul du vecteur directeur unitaire de l'axe (U2) GOTO 1000 ENDIF YL=1.D0/SQRT(YL) DO i=1,IDIM U2(i)=U2(i)*YL ENDDO ENDIF GOTO 449 C On a trouve le mot DIREction 44804 IDIREC=1 IF (IRETOU.EQ.0) THEN * WRITE(IOIMP,*) 'Lecture dun chpoint' IF (IERR.NE.0) GOTO 1000 IDICHP=1 ELSE j=(KPOINT-1)*idimp1 YL=0.D0 DO i=1,IDIM XNOR(i)=XCOOR(j+i) YL=YL+XNOR(i)*XNOR(i) ENDDO GOTO 1000 ENDIF YL=1.D0/SQRT(YL) DO i=1,IDIM XNOR(i)=XNOR(i)*YL ENDDO ENDIF GOTO 449 C Lecture eventuelle de DDLs : C ------------------------------ 445 CONTINUE C La liste des grandeurs PRIMALEs/DDLs est un LISTMOTS. cbp,2020-11-12 : possibilite non documentee => ajout a la notice c + ajout de qq lignes pour rendre listmot dual facultatif IF (IERR.NE.0) GOTO 1000 IF (IRETOU.EQ.0) GOTO 446 IF (IERR.NE.0) GOTO 1000 C -cas ou la liste des grandeurs DUALEs n'est pas fournie. IF(IRETO2.EQ.0) THEN SEGACT,MLMOT1 IF (IMOT1.LE.0) THEN SEGDES,MLMOT1 GOTO 1000 ENDIF DO IMOT = 1, IMOT1 IF (IERR.NE.0.or.idd.LE.0) THEN GOTO 1000 ENDIF MOTDDL(**)=NOMDU(idd) ENDDO SEGDES,MLMOT1 C -cas ou la liste des grandeurs DUALEs est un LISTMOTS fourni. cbp,2020-11-12 : on ne documente pas cette possibilite car pas le temps ELSE SEGACT,MLMOT1,MLMOT2 IF (IMOT1.LE.0) THEN SEGDES,MLMOT1,MLMOT2 GOTO 1000 SEGDES,MLMOT1,MLMOT2 GOTO 1000 ENDIF DO IMOT = 1, IMOT1 ENDDO SEGDES,MLMOT1,MLMOT2 ENDIF GOTO 449 446 CONTINUE C La liste des grandeurs PRIMALEs/DDLs sont des MOTs. C La liste des ddls autorises = NOMDD defini dans BDATA.ESO C On doit lire au moins 1 ddl (car sinon MOTDDL est vide !) IF (IERR.NE.0) GOTO 1000 IF (IMOT.NE.0) THEN MOTDDL(**)=NOMDD(IMOT) MOTDDL(**)=NOMDU(IMOT) ELSE IF (IERR.NE.0) GOTO 1000 IF (IMOT.EQ.0) GOTO 447 MOTDDL(**)=CHADDL MOTDDL(**)=' ' ENDIF GOTO 446 C Conditions aux limites via un MODELE --> non documente, non plus ! 447 CONTINUE IBDDL=MOTDDL(/2) IPMODL = 0 IF (IERR.NE.0) GOTO 1000 IF (IPMODL.EQ.0) THEN DO i = 2, IBDDL, 2 IF (MOTDDL(i).EQ.' ') THEN MOTERR(1:4)=MOTDDL(i-1) ENDIF ENDDO ELSE IF (IERR.NE.0) GOTO 1000 IF (IERR.NE.0) GOTO 1000 SEGACT,MLMOT1,MLMOT2 DO i = 2, IBDDL, 2 IF (MOTDDL(i).EQ.' ') THEN IF (iplac.EQ.0) THEN MOTERR(1:4)=MOTDDL(i-1) ELSE ENDIF ENDIF ENDDO SEGSUP,MLMOT1,MLMOT2 ENDIF IF (IERR.NE.0) GOTO 1000 449 IBDDL=MOTDDL(/2) C Verification que le nombre de DDLs a bloquer n'est pas nul C IF (IBDDL.EQ.0) GOTO 1000 C----------------------------------------------------------------------- C Fin de la lecture des mots (DEPL,ROTA...) ou des DDLs C----------------------------------------------------------------------- C Recherche du maillage MELEME de type POINT : C ---------------------------------------------- KOBJET=0 MELEME=0 C On cherche d'abord si on a un POINT que l'on transformera en POI1 C sinon on cherche un maillage que l'on transforme en POI1 si besoin IF (IERR.NE.0) GOTO 1000 IF (IRETOU.EQ.0) THEN IF (IERR.NE.0) GOTO 1000 MELEME=KOBJET SEGACT,MELEME NBPOIN=NUM(/2) ELSE MELEME=KPOINT NBPOIN=1 ENDIF * IF (IDICHP.EQ.1) THEN * Cas ou la direction est donnee par un chpoint * On construit un segment MTRAV contenant les directions normees * WRITE(IOIMP,*) 'Transfo du chpoint en MTRA' * 1) Recopie des composantes et valeurs pertinentes du chpoint * dans le TMTRAV * Composantes recherchees (LMDEPL.MLMOTS) * sur les points du maillage MELEME JGN=4 JGM=IDIM SEGINI,LMDEPL DO I=1,IDIM ENDDO SEGSUP LMDEPL SEGACT MELEME IF (IRET.NE.0) THEN GOTO 1000 ENDIF * 2) Calcul des normes * WRITE(IOIMP,*) 'Calcul des normes ' SEGACT MTRAV*MOD DO IBPOIN=1,NBPOIN YL=0.D0 DO I=1,IDIM XNOR(I)=BB(I,IBPOIN) YL=YL+XNOR(I)*XNOR(I) ENDDO IF (YL.LT.XPETIT) THEN GOTO 1000 ENDIF YL=1.D0/SQRT(YL) DO I=1,IDIM BB(I,IBPOIN)=XNOR(I)*YL ENDDO ENDDO * WRITE(IOIMP,*) 'Fin du calcul des normes ' * SEGPRT,MTRAV ENDIF C Determination du nombre de multiplicateurs NNMAT par noeud de MELEME C NNMAT correspond au nombre de DDLs a bloquer par noeud = nombre de C multiplicateurs a creer par noeud (1 multiplicateur) = NRIGEL C Dans les cas RADIal et DIREction, on a une seule matrice par noeud. C Dans les autres cas, autant de matrices que MOTDDL(/2)/2. NNMAT=1 IF (IDIREC+IRADIA.EQ.0) NNMAT=IBDDL/2 C Initialisation de l'objet RIGIDITE associe aux BLOCAGES NRIGE=8 NRIGEL=NNMAT SEGINI,MRIGID MTYMAT=MOTRIG IFORIG=IFOUR ICHOLE=0 IMGEO1=0 IMGEO2=0 KRIGI=MRIGID C A chaque multplicateur est associe un nouveau noeud C NBPOIN : nombre de points du maillage MELEME a bloquer NBNO=nbpts NBNOI=NBNO SEGADJ,MCOORD C Boucle sur le nombre de DDLs a bloquer DO IAA=1,NNMAT C Pour chaque noeud du maillage, creation d'un noeud associe (place C initialement a l'origine) au IAA-eme DDL DO i=1,NBPOIN DO j=1,idimp1 XCOOR(NBNOI*idimp1+j)=0.D0 ENDDO NBNOI=NBNOI+1 ENDDO C Creation du maillage MELEME de MULTiplicateurs associe aux BLOCAGES NBSOUS=0 NBREF=0 NBNN=2 NBELEM=NBPOIN SEGINI,IPT1 KIPT1=IPT1 IPT1.ITYPEL=22 DO i=1,NBPOIN j=(IAA-1)*NBPOIN+i IPT1.NUM(2,i)=NUM(1,i) IPT1.ICOLOR(i)=IDCOUL C Correction : Les pts mult sont a la meme position que les noeuds IREF3=(NUM(1,i)-1)*idimp1 DO j=1,IDIM XCOOR(IREF1+j)=XCOOR(IREF3+j) ENDDO ENDDO C Creation des RAIDEURS associees au IAA-eme multplicateur (DDL) IRIGEL(1,IAA)=KIPT1 IRIGEL(2,IAA)=0 IRIGEL(5,IAA)=NIFOUR IRIGEL(6,IAA)=NILATE NLIGRP=2 IF (IDIREC+IRADIA.NE.0) NLIGRP=1+LDEPL NLIGRD=NLIGRP C Remplissage du tableau des DESCripteurs de RIG SEGINI,DESCR IRIGEL(3,IAA)=DESCR NOELEP(1)=1 NOELEP(2)=2 NOELED(1)=1 NOELED(2)=2 IF (IDIREC+IRADIA.EQ.0) THEN j=2*(IAA-1) LISINC(1)='LX ' LISINC(2)=MOTDDL(j+1) LISDUA(1)='FLX ' LISDUA(2)=MOTDDL(j+2) ELSE DO i=1,LDEPL NOELEP(1+i)=2 NOELED(1+i)=2 IF (IROTA.NE.1) THEN LISINC(1+i)=MODEPL(IADEPL+i) LISDUA(1+i)=MODEDU(IADEPL+i) ELSE LISINC(1+i)=MOROTA(IADEPL+i) LISDUA(1+i)=MORODU(IADEPL+i) ENDIF ENDDO LISINC(1)='LX ' LISDUA(1)='FLX ' ENDIF SEGDES DESCR C** NLIGRP=3 C** IF (IDIREC+IRADIA.NE.0) NLIGRP=IDIM+2 C** NLIGRD=NLIGRP NELRIG=NBPOIN SEGINI,XMATRI IRIGEL(4,IAA)=XMATRI COERIG(IAA)=1.D0 C Remplissage de la matrice de rigidite RE : C Il faut distinguer les cas IRADIA, IDIREC et les autres C IRADIA : Il faut calculer la DIREction puis identique a IDIREC C IDIREC : La DIRECTION est stockee dans le vecteur norme XNOR C AUTRES : La matrice est predefinie dans RIG C Option RADIAL : Calcul prealable de la direction pour chaque noeud IF (IRADIA.NE.0) THEN DO IB=1,NBPOIN j=(NUM(1,IB)-1)*idimp1 DO i=1,IDIM XNOR(i)=XCOOR(j+i)-U1(i) ENDDO IF (IDIM.EQ.2) THEN YL=XNOR(1)*XNOR(1)+XNOR(2)*XNOR(2) RETURN ENDIF YL=1.D0/SQRT(YL) IF (IRADIA.EQ.1) THEN XNOR(1)=XNOR(1)*YL XNOR(2)=XNOR(2)*YL ELSE IF (IRADIA.EQ.2) THEN XX=XNOR(1) XNOR(1)=-XNOR(2)*YL XNOR(2)=XX*YL ENDIF ELSE YL=XNOR(1)*U2(1)+XNOR(2)*U2(2)+XNOR(3)*U2(3) XL=0.D0 DO i=1,3 XNOR(i)=XNOR(i)-YL*U2(i) XL=XL+XNOR(i)*XNOR(i) ENDDO RETURN ENDIF IF (IRADIA.EQ.1) THEN XL=1.D0/SQRT(XL) XNOR(1)=XNOR(1)*XL XNOR(2)=XNOR(2)*XL XNOR(3)=XNOR(3)*XL ELSE IF (IRADIA.EQ.2) THEN XX=XNOR(1) YY=XNOR(2) ZZ=XNOR(3) XNOR(1)=YY*U2(3)-ZZ*U2(2) XNOR(2)=ZZ*U2(1)-XX*U2(3) XNOR(3)=XX*U2(2)-YY*U2(1) ENDIF ENDIF C XNOR contient la direction normee RE(1,1,IB)=0.D0 RE(2,1,IB)=XNOR(1) RE(3,1,IB)=XNOR(2) RE(1,2,IB)=RE(2,1,IB) RE(1,3,IB)=RE(3,1,IB) IF (IDIM.EQ.3) THEN RE(4,1,IB)=XNOR(3) RE(1,4,IB)=RE(4,1,IB) ENDIF ENDDO C Option DIRECTION ELSE IF (IDIREC.EQ.1) THEN IF (IDICHP.EQ.1) THEN DO I=1,NBPOIN RE(1,1,i)=0.D0 RE(2,1,i)=BB(1,I) RE(3,1,i)=BB(2,I) RE(1,2,i)=RE(2,1,i) RE(1,3,i)=RE(3,1,i) RE(2,3,i)=RE(3,2,i) IF (IDIM.EQ.3) THEN RE(4,1,i)=BB(3,I) RE(1,4,i)=RE(4,1,i) ENDIF ENDDO ELSE DO i=1,NBPOIN RE(1,1,i)=0.D0 RE(2,1,i)=XNOR(1) RE(3,1,i)=XNOR(2) RE(1,2,i)=RE(2,1,i) RE(1,3,i)=RE(3,1,i) RE(2,3,i)=RE(3,2,i) IF (IDIM.EQ.3) THEN RE(4,1,i)=XNOR(3) RE(1,4,i)=RE(4,1,i) ENDIF * DO i=1,NBPOIN ENDDO SEGDES,xMATRI ENDIF C Autres options : ELSE DO i=1,NBPOIN RE(1,1,i)=0.D0 RE(2,1,i)=1.D0 RE(2,2,i)=0.D0 RE(1,2,i)=RE(2,1,i) * DO i=1,NBPOIN ENDDO SEGDES,xMATRI SEGDES,IPT1 ENDIF SEGDES,xMATRI ENDDO C Fin de la boucle sur les IAA DDLs a bloquer IF (IDICHP.EQ.1) THEN * WRITE(IOIMP,*) 'Destruction du segment de travail' SEGSUP,MTRAV ENDIF C Fin normale du traitement krigi=mrigid SEGDES,MRIGID IF (MELEME.GT.0) THEN IF (KOBJET.EQ.0) THEN SEGSUP,MELEME ELSE IF (MELEME.EQ.KOBJET) THEN SEGDES,MELEME ELSE ** SEGSUP,MELEME SEGDES,MELEME ENDIF ENDIF ENDIF C Il peut rester des maillages actifs MRIGID=KRIGI SEGACT,MRIGID DO i=1,IRIGEL(/2) MELEME=IRIGEL(1,i) IF (MELEME.GT.0) SEGDES,MELEME ENDDO SEGDES,MRIGID C Sortie du sousprogramme (normale ou en cas d'erreur) 1000 SEGSUP,MSWBLO SEGDES,MCOORD END
© Cast3M 2003 - Tous droits réservés.
Mentions légales