C RELSUR SOURCE PV090527 26/04/30 21:16:11 12529 C RELSUR SOURCE BP208322 17/04/18 21:15:13 9395 SUBROUTINE RELSUR(IMODEL, MRIGID) C*********************************************************************** C cet operateur crée une matrice élémentaire de rigidité c pour imposer les relation portées par un modele de sure C C ENTREES : C ________ C C IMODEL pointeur sur le modele élémentaire C C ENTREES/SORTIE : C ________ C c MRIGID rigidité chapeau dans laquelle on va écrire c (à la suite des sous matrices déja présentes) les rigidités c voulues. C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C -INC PPARAM -INC CCOPTIO -INC SMRIGID C-INC SMINTE -INC SMMODEL -INC SMELEME -INC SMCHAML -INC CCHAMP -INC CCGEOME -INC SMCOORD C C C Petit tableau des "couleurs" des relations de conformite DIMENSION LCOLOR(6) DATA LCOLOR / 1, 3, 6, 10, 16, 24 / C C nombre de sous matrice de mrigid (va être ammené a changé) NRIGEL = MRIGID.IRIGEL(/2) IPT1=IMODEL.imamod SEGACT, IPT1 IDEBUT = LCOLOR(IPT1.ICOLOR(1)) - 3 c récupérations du nom des composantes obligatoire du modele de sure nomid=IMODEL.lnomid(1) SEGACT, nomid SEGACT,MCOORD*MOD C********************************************************************** C Boucle sur les composantes primales obligatoires du sure C********************************************************************** ICOMP=0 DO 2 ICOMP=1,nomid.lesobl(/2) C==================== c creation d'un maillage avec un premier noeud par lélément c correspondant à un multiplicateur de lagrange C==================== NBNN=IPT1.NUM(/1)+1 NBELEM=IPT1.NUM(/2) NBSOUS=0 NBREF=0 SEGINI, IPT4 IPT4.ITYPEL=22 DO 1 J=1, NBELEM ipt4.icolor(j)=IPT1.icolor(j) DO 11 I=1,IPT1.NUM(/1) IPT4.NUM(I+1,J)=IPT1.NUM(I,J) 11 CONTINUE 1 CONTINUE C creation n'un nouveau noeud pour supporter chaque multiplicateur de lagrange NBPT1= nbpts NBPTS=NBPT1+(IPT4.NUM(/2)) SEGADJ MCOORD DO 12 J=1, NBELEM NGLOB=(NBPT1+J-1)*(IDIM+1) C remplissage des coordonees des nouveux points DO 13 ID= 1,IDIM XCOOR(NGLOB+ID)=XCOOR((IPT4.NUM(2,J)-1)*(IDIM+1)+ID) 13 CONTINUE IPT4.NUM(1,J) = NBPT1 + J 12 CONTINUE C==================== C *** SEGMENT XMATRI C==================== NLIGRD=IPT4.NUM(/1) NLIGRP=NLIGRD NELRIG=IPT4.NUM(/2) rigrel=0 SEGINI, XMATRI DO 3 I=1,NELRIG RE(1,2,i)=-1. RE(2,1,i)=-1. DO 4 ICAZ=3,NLIGRD RE(1,ICAZ,i)=XCOEFF(IDEBUT+ICAZ) RE(ICAZ,1,i)=RE(1,ICAZ,i) 4 CONTINUE 3 CONTINUE C write(*,*) 'COMPOSANTE OBLIGATOIRE DU SURE :' C write(*,*) (nomid.lesobl(ICOMP)) C write(*,*) 'MATRICE elementaire :' C DO 5 I=1,NLIGRD C write(*,*) (RE(i,iou,1), iou=1,NLIGRD) C 5 CONTINUE C==================== C *** SEGMENT DESCR C==================== NEXIST=0 DO 6 I=1,LNOMDD IF (NOMDD(I).EQ.nomid.lesobl(ICOMP)) NEXIST = I 6 CONTINUE IF (NEXIST.EQ.0) THEN CALL ERREUR(837) ENDIF SEGINI, DESCR LISINC(1)='LX ' LISDUA(1)='FLX ' NOELEP(1)=1 NOELED(1)=1 DO 7 I=2,NLIGRD LISINC(I)=NOMDD(NEXIST) LISDUA(I)=NOMDU(NEXIST) NOELEP(i)=i NOELED(i)=i 7 CONTINUE C==================== C *** stockage de la rigidité construite dans MRIGID C==================== NRIGEL = NRIGEL+1 SEGADJ, MRIGID SEGACT, MRIGID*mod COERIG(NRIGEL)=1. IRIGEL(1,NRIGEL)=IPT4 IRIGEL(3,NRIGEL)=DESCR IRIGEL(4,NRIGEL)=XMATRI c SEGDES, MRIGID SEGDES, DESCR SEGDES, XMATRI SEGDES, IPT4 C********************************************************************** C FIN Boucle sur les composantes primales obligatoires du sure C********************************************************************** 2 CONTINUE c write (*,*) 'nomid.lesfac(/2)', nomid.lesfac(/2) IF (nomid.lesfac(/2).EQ.0) goto 100 c récupération du champ d'enrichissement c MCHEX1= IMODEL.IVAMOD(1) c SEGACT, MCHEX1 c MCHAM1= MCHEX1.ICHAML(1) MCHAM1= IMODEL.IVAMOD(1) SEGACT, MCHAM1 C********************************************************************** C Boucle sur les composantes primales facultatives du sure C********************************************************************** ICOMP=0 DO 101 ICOMP=1,nomid.lesfac(/2) C++++ choix du type d'enrichisement de la composante ICOMP IF (nomid.lesfac(ICOMP).EQ.'AX') MELVA1=MCHAM1.IELVAL(1) IF (nomid.lesfac(ICOMP).EQ.'AY') MELVA1=MCHAM1.IELVAL(1) IF (nomid.lesfac(ICOMP).EQ.'AZ') MELVA1=MCHAM1.IELVAL(1) IF (nomid.lesfac(ICOMP).EQ.'B1X') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'B1Y') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'B1Z') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'C1X') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'C1Y') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'C1Z') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'D1X') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'D1Y') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'D1Z') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'E1X') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'E1Y') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'E1Z') MELVA1=MCHAM1.IELVAL(2) IF (nomid.lesfac(ICOMP).EQ.'B2X') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'B2Y') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'B2Z') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'C2X') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'C2Y') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'C2Z') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'D2X') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'D2Y') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'D2Z') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'E2X') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'E2Y') MELVA1=MCHAM1.IELVAL(3) IF (nomid.lesfac(ICOMP).EQ.'E2Z') MELVA1=MCHAM1.IELVAL(3) SEGACT MELVA1 C==================== c creation d'un maillage de multiplicateurs de lagranges enrichis C==================== NBNN=IPT1.NUM(/1)+1 NBELEM=IPT1.NUM(/2) NBSOUS=0 NBREF=0 SEGINI, IPT4 IPT4.ITYPEL=22 IELENR=0 c++++ BOUCLE sur les éléments de ipt1 DO 102 JELEM=1,IPT1.NUM(/2) NEXIST=0 ipt4.icolor(jelem)=IPT1.icolor(jelem) JNUM = IPT1.NUM(1,JELEM) c+++ Recherche d'une valeur non nulle du champ d'enrichissement VENR1 = MELVA1.VELCHE(1,JELEM) C On prend les elements dont le hanging node est enrichi IF (VENR1.GT.0) THEN NEXIST=NEXIST+1 C On prend les element dont tout les autres noeuds sont enrichis ELSE DO 121 JNOEUD= 2 , IPT1.NUM(/1) VENR1 = MELVA1.VELCHE(JNOEUD,JELEM) IF (VENR1.eq.0) GOTO 121 NEXIST=NEXIST+1 121 CONTINUE ENDIF IF (nexist.eq.0) GOTO 102 IELENR= IELENR+1 C On recopie dans IPT4 les elements de ipt1 sur lequel on veux c imposer une relation de compatibilité DO 122 I=1,IPT1.NUM(/1) IPT4.NUM(I+1,IELENR)=IPT1.NUM(I,JELEM) 122 CONTINUE 102 CONTINUE NBELEM=IELENR SEGADJ IPT4 IF (ielenr.eq.0) then segsup ipt4 goto 101 endif C======================================================================= C creation n'un nouveau noeud pour supporter chaque multiplicateur de lagrange C======================================================================= NBPT1=nbpts NBPTS=NBPT1+IELENR SEGADJ,MCOORD DO 103 J=1,IPT4.NUM(/2) NGLOB=(NBPT1+J-1)*(IDIM+1) C remplissage des coordonees des nouveux points DO 131 ID= 1,IDIM XCOOR(NGLOB+ID)=XCOOR((IPT4.NUM(2,J)-1)*(IDIM+1)+ID) 131 CONTINUE IPT4.NUM(1,J) = NBPT1 + J 103 CONTINUE C==================== C *** SEGMENT XMATRI C==================== NLIGRD=IPT4.NUM(/1) NLIGRP=NLIGRD NELRIG=IPT4.NUM(/2) rigrel=0 SEGINI, XMATRI c++++ BOUCLE sur les éléments de ipt4 DO 104 IELEM=1,NELRIG RE(1,2,IELEM)=-1. RE(2,1,IELEM)=-1. DO 141 ICAZ=3,NLIGRD RE(1,ICAZ,IELEM)=XCOEFF(IDEBUT+ICAZ) RE(ICAZ,1,IELEM)=RE(1,ICAZ,IELEM) 141 CONTINUE 104 CONTINUE C write(*,*) 'COMPOSANTE FACULTATIVE DU SURE :' C write(*,*) (nomid.lesfac(ICOMP)) C write(*,*) 'MATRICE elementaire :' C DO 6666 I=1,NLIGRD C write(*,*) (RE(i,iou,1), iou=1,NLIGRD) C 6666 CONTINUE C==================== C *** SEGMENT DESCR C==================== NEXIST=0 DO 105 ICO1=1,LNOMDD IF (NOMDD(ICO1).EQ.nomid.lesfac(ICOMP)) NEXIST = ICO1 105 CONTINUE IF (NEXIST.EQ.0) THEN CALL ERREUR(837) ENDIF SEGINI, DESCR LISINC(1)='LX ' LISDUA(1)='FLX ' NOELEP(1)=1 NOELED(1)=1 DO 106 ICO2=2,NLIGRD LISINC(ICO2)=NOMDD(NEXIST) LISDUA(ICO2)=NOMDU(NEXIST) NOELEP(ico2)=ico2 NOELED(ico2)=ico2 106 CONTINUE C==================== C *** stockage de la rigidité construite dans MRIGID C==================== C Ajustement du segment rigidite NRIGEL=NRIGEL +1 SEGADJ, MRIGID C* ICHOLE=0 C* IMGEO1=0 C* IMGEO2=0 C* IFORIG=0 C* ISUPEQ=0 COERIG(NRIGEL)=1. IRIGEL(1,NRIGEL)=IPT4 IRIGEL(3,NRIGEL)=DESCR IRIGEL(4,NRIGEL)=XMATRI SEGDES, DESCR SEGDES, XMATRI SEGDES, IPT4 C********************************************************************** C FIN Boucle sur les composantes primales facultatives du sure C********************************************************************** 101 CONTINUE 100 CONTINUE SEGDES, nomid segdes,IPT1, MCOORD C RETURN END