C RESOU SOURCE PV090527 23/03/15 21:15:02 11630 SUBROUTINE RESOU IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C C **** CET OPERATEUR SERT A RESOUDRE UN SYSTEME D EQUATIONS LINEAIRES C **** CHPOINT = RESOU RIGIDITE CHPOINT C C -INC PPARAM -INC CCOPTIO -INC SMRIGID -INC SMTEXTE -INC SMTABLE -INC SMCHPOI -INC SMELEME -INC SMLCHPO -INC CCREEL PARAMETER(ZERO=0.D0) SEGMENT IDEMEM(0) segment ideme0(idemem(/1),30) segment ideme1(idemem(/1),30) segment idnote(0) C CHARACTER*4 LISM(9) CHARACTER*(8) CHARRE1 CHARACTER*72 CHARRE REAL*8 XVA LOGICAL ILOG,ILUG,casfimp DATA LISM/'NOID','NOUN','ENSE','GRAD','CHOL','STAB','ELIM', >'NOST','SOUC'/ DATA ILOG/.FALSE./ C C------------------------------------------------------- c LECTURE ET INITIALISATION c LECTURE DES OPTIONS XVA=REAL(0.D0) IOB=0 * sauvegarde norinc * write(6,*) 'norinc norind ',norinc,norind norins=norinc iverif=1 ipt8=0 iunil=0 * le defaut est de faire une passe d'elimination nelim=30 * experimentalement 2 passes est mieux nelim=2 IMTVID=0 NOUNIL=0 NOID=0 NOEN=1 IGRADJ = 0 ICHSKI = 0 INSYM = 0 KIKI=0 KSYMET = 0 IPSHPO = 0 ISTAB=0 ISOUCI = 0 ifochs = -99 5 CONTINUE CALL LIRMOT(LISM,9,KIKI,0) IF (KIKI.EQ.1) NOID=1 IF (KIKI.EQ.2) NOUNIL=1 IF (KIKI.EQ.3) NOEN=0 * IF (KIKI.EQ.4) IGRADJ = 1 * IF (KIKI.EQ.5) ICHSKI = 1 * IF (KIKI.EQ.4.OR.KIKI.EQ.5) KSYMET = 1 IF (KIKI.eQ.6) ISTAB=1 IF (KIKI.eQ.7) then call lirent(nelim,1,iretou) nelim=min(30,max(0,nelim)) endif IF (KIKI.eQ.8) ISTAB=0 IF (KIKI.eQ.9) ISOUCI=1 IF (KIKI.NE.0) GOTO 5 if (noid.eq.1) iverif=0 IF(NUCROU.EQ.0) THEN ICHSKI=1 ELSEIF(NUCROU.EQ.1) THEN IGRADJ=1 KSYMET=1 ENDIF * WRITE(6,*) ' nucrou', nucrou * IF ( IGRADJ + ICHSKI .EQ. 0 ) ICHSKI = 1 c LECTURE DE LA RIGIDITE CALL LIROBJ('RIGIDITE',IPOIRI,1,IRETOU) IF(IERR.NE.0) GO TO 5000 IPRIGO=IPOIRI C c LECTURE DE LA PRECISION PREC=REAL(xspeti) CALL LIRREE(PREC,0,IRETOU) IF(IERR.NE.0) GO TO 5000 C REMPLISSAGE DU 2ND MEMBRE IDEMEM(**) A PARTIR DE ... c ... 'CHPOINT' SEGINI IDEMEM 1 CONTINUE CALL LIROBJ('CHPOINT ',ISECO,0,IRETOU) IF(IRETOU.NE.0) THEN IDEMEM(**)=ISECO mchpoi=iseco segact mchpoi * write(6,*) ' extension idemem 1 ',idemem(/1) GO TO 1 ENDIF c ... 'TABLE DE SOUS-TYPE LIAISONS_STATIQUES' CALL LIRTAB('LIAISONS_STATIQUES',ITBAS,0,IRET) c ... 'LISTCHPO' CALL LIROBJ('LISTCHPO',ISECO,0,IRETOU) IF(IRETOU.NE.0) THEN mlchpo=ISECO segact mlchpo n1 = ichpoi(/1) do iu = 1 , n1 idemem(**) = ichpoi(iu) * write(6,*) ' extension idemem 2 ',idemem(/1) enddo segdes mlchpo segini mlchpo ipshpo = mlchpo ENDIF IF (IERR.NE.0) RETURN IF (ITBAS.NE.0 .AND. IIMPI.EQ.333) THEN WRITE(IOIMP,*) 'on a lu la table des conditions aux limites' ENDIF if (itbas.ne.0) then mtab1 = itbas segact mtab1 ima = mtab1.mlotab - 1 segini idnote im = 0 segdes mtab1 else goto 90 endif * boucle en cas de résolutions successives avec table 80 continue im = im + 1 itmod = 0 ichp0 = 0 if (im.gt.ima) then if (idemem(/1).gt.0) goto 90 * pas de champs de force call erreur(1) return endif CALL ACCTAB(ITBAS,'ENTIER',IM,0.d0,' ',.true.,IP0, & 'TABLE',I1,X1,CHARRE,.true.,ITMOD) if (ierr.ne.0) return c table itmod trouvee --> on recupere la force if (itmod.gt.0) then CALL ACCTAB(ITMOD,'MOT',0,0.d0,'FORCE',.true.,IP0, & 'CHPOINT',I1,X1,CHARRE,.true.,ICHP0) if (ierr.ne.0) return if (ichp0.gt.0) then idemem(**) = ichp0 * write(6,*) ' extension idemem 3 ',idemem(/1) idnote(**) = im else call erreur(1) return endif c on cree le point repere ici CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN) CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0, & 'POINT',0,0.0D0,' ',.TRUE.,IPOIN) endif goto 80 IF (IERR.NE.0) RETURN C------------------------------------------------------- c DEBUT DU TRAVAIL 90 continue segini ideme0,ideme1 * verification pas de blocage en double *** call verlag(ipoiri) if (ierr.ne.0) return * y a t il des matrices de relations non unilaterales ipoir0 = ipoiri mrigid=ipoiri C call prrigi(ipoiri,1) segact mrigid ifochs = iforig nrige= irigel(/1) idepe=0 * write(ioimp,*) 'dans resou mrigid iforig ',mrigid,iforig nbr = irigel(/2) do 1000 irig = 1,nbr meleme=irigel(1,irig) segact meleme if ((irigel(6,irig).eq.0.or.nounil.eq.1).and.itypel.eq.22) > idepe=idepe+num(/2) if (irigel(6,irig).ne.0) iunil=1 if (irigel(6,irig).eq.2) nelim=30 if(jrcond.ne.0) nelim=30 if (irigel(7,irig).ne.0) then insym=1 ichski=0 endif 1000 continue * elimination recursive des conditions aux limites * on la fait en gradient conjugue ou en appel de unilater nfois=nelim-1 if (igradj.eq.1.or.(iunil.eq.1.and.nounil.eq.0)) nfois=29 lagdua=0 imult=1 icond=idepe icondi=(icond*10)/9+1 if=0 do ifois=1,nfois if(imult.ne.0.and.icond.ne.0.and.(icond*10)/9.lt.icondi.and. > (icondi-icond.gt.0.or.igradj.eq.1)) then icondi=icond lagdua=-1 if=if+1 if(ierr.ne.0) return call resouc(mrigid,mrigic,idemem,ideme0,ideme1, > nounil,lagdua,icond,imult,if,imtvid,nelim) ** write(ioimp,*) ' passe ',if,' condition ',icond,ifois if(ierr.ne.0) return mrigid=mrigic endif enddo * Si on n'a pas reussi a tout eliminer, on fait encore une passe pour creer lagdua lagdua=0 if (iunil.eq.0.or.nounil.eq.1) then if (icond.ne.0) then if=if+1 if(ierr.ne.0) return call resouc(mrigid,mrigic,idemem,ideme0,ideme1, > nounil,lagdua,icond,imult,if,imtvid,nelim) * write(ioimp,*) ' passe ','finale',' condition ',icond if(ierr.ne.0) return mrigid=mrigic endif endif ** write (ioimp,*) 'nombre de passes',if if (idepe.ne.0) noid = 1 ipoiri=mrigid * call prrigi(ipoiri,1) C------------------------------------------------------- * * Si au moins une des matrices n'est pas symétrique, on passera * par le solveur non-symétrique LDMT. * SEGACT MRIGID*MOD NRG = IRIGEL(/1) NBR = IRIGEL(/2) C ... Ceci peut arriver si par exemple on extrait la partie C symétrique d'une matrice purement antisymétrique ... * IF(NBR.EQ.0) THEN * SEGDES MRIGID * CALL ERREUR(727) * RETURN * ENDIF C ... Mais avant on va tester si la normalisation des variables C primales et duales a été demandée - ceci entraîne la perte C de la symétrie ... IF(NORINC.GT.0 .AND. NORIND.GT.0) THEN IF(KSYMET.EQ.1) THEN CALL ERREUR(19) SEGDES,MRIGID RETURN ENDIF INSYM = 1 IGRADJ = 0 ICHSKI = 0 GOTO 15 ENDIF IF (NRG.GE.7) THEN C ... On teste si la matrice contient des matrices non symétriques ... C ... Si OUI, ce n'est pas la peine de faire les autres tests ... DO 9 IN = 1,NBR IANTI=IRIGEL(7,IN) IF(IANTI.GT.0) THEN C ... On vérifie si l'utilisateur n'a pas demandé explicitement C la résolution par Choleski ou gradient conjugué, C si OUI on râle puis on s'en va !!! ... IF(KSYMET.EQ.1) THEN CALL ERREUR(1126) SEGDES,MRIGID RETURN ENDIF IF(NORINC.NE.0.AND.NORIND.EQ.0) THEN * on supprime la normalisation au lieu de faire une erreur norinc=0 ** CALL ERREUR(760) ** SEGDES,MRIGID ** RETURN ENDIF INSYM = 1 IGRADJ = 0 ICHSKI = 0 GOTO 15 ENDIF 9 CONTINUE ENDIF 15 CONTINUE C C **** ON S'ASSURE QU'IL N'Y A PAS D'APPUIS UNILATERAUX C if (iunil.eq.0) goto 30 IF(IRIGEL(/1).LT.6) GO TO 30 IF (NOUNIL.EQ.1) GOTO 30 21 CONTINUE C C **** EXISTENCE DES APPUIS UNILATERAUX C **** SI ON EST DEJA PASSE DANS LES APPUIS UNILATERAUX C ISUPEQ POINTE SUR UNE TABLE CONTENANT LES DONNEES A PASSER C A LA PROCEDURE UNILATER C ISUPLO=ISUPEQ IF (ISUPLO.NE.0) GOTO 27 NNOR=0 DO 22 I=1,IRIGEL(/2) IF(IRIGEL(6,I).EQ.0) NNOR=NNOR+1 22 CONTINUE IF(NNOR.EQ.0) THEN CALL ERREUR(312) norinc=norins RETURN ENDIF NRIGE=IRIGEL(/1) NRIGEL=NNOR SEGINI RI1 RI1.IFORIG = IFORIG c* RI1.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe RI1.MTYMAT = ' ' NRIGEL=IRIGEL(/2)-NNOR SEGINI RI2 RI2.IFORIG = IFORIG c* RI2.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe RI2.MTYMAT = ' ' II1=0 II2=0 DO 23 I=1,IRIGEL(/2) IF(IRIGEL(6,I).NE.0) THEN RI3=RI2 II2=II2+1 II=II2 ELSE RI3=RI1 II1=II1+1 II=II1 ENDIF DO 24 J=1,NRIGE RI3.IRIGEL(J,II) = IRIGEL(J,I) 24 CONTINUE RI3.COERIG(II)=COERIG(I) 23 CONTINUE * RI1 raideur sans condition unilaterale * RI2 conditions unilaterales CALL CRTABL(MTABLE) ISUPEQ=MTABLE * il faut aussi mettre isupeq dans la raideur initiale if (jrsup.ne.0) mrigid=jrsup segact mrigid iri1s=jrelim iri2s=mrigid MRIGID=IPRIGO SEGACT MRIGID*MOD ISUPEQ=MTABLE if (idepe.ne.0) then * il faut extraire de la matrice initiale (ipoir0) les conditions unilaterales mrigid=iri2s segact mrigid nrigel=0 do 40 i=1,irigel(/2) if (irigel(6,i).eq.0) nrigel=nrigel+1 40 continue if (nrigel.eq.0) call erreur(312) if (ierr.ne.0) then norinc=norins return endif nrige=irigel(/1) segini ri4 ri4.iforig = iforig ri4.mtymat = mtymat ii1=0 nrigel=irigel(/2)-nrigel segini ri5 ri5.iforig = iforig ri5.mtymat = mtymat ii2=0 do 41 i=1,irigel(/2) if (irigel(6,i).ne.0) goto 42 ii1=ii1+1 do j=1,nrige ri4.irigel(j,ii1)=irigel(j,i) enddo ri4.coerig(ii1)=coerig(i) goto 41 42 continue ii2=ii2+1 do j=1,nrige ri5.irigel(j,ii2)=irigel(j,i) enddo ri5.coerig(ii2)=coerig(i) 41 continue segdes mrigid,ri4 endif ri3=iri1s * segact ri1,ri2,ri3,ri4,ri5 CALL ECCTAB(MTABLE,'ENTIER ',1,XVA,' ',ILOG,IOB, $ 'RIGIDITE',IOB,XVA,' ',ILOG,RI1) CALL ECCTAB(MTABLE,'ENTIER ',2,XVA,' ',ILOG,IOB, $ 'RIGIDITE',IOB,XVA,' ',ILOG,RI2) CALL ECCTAB(MTABLE,'ENTIER ',3,XVA,' ',ILOG,IOB, $ 'LOGIQUE ',IOB,XVA,' ',ILOG,IOB) ** if(idepe.ne.0) then ** CALL ECCTAB(MTABLE,'ENTIER ',8,XVA,' ',ILOG,IOB, ** $ 'RIGIDITE',IOB,XVA,' ',ILOG,iri1s) ** CALL ECCTAB(MTABLE,'ENTIER ',9,XVA,' ',ILOG,IOB, ** $ 'RIGIDITE',IOB,XVA,' ',ILOG,ri4 ) ** CALL ECCTAB(MTABLE,'ENTIER ',12,XVA,' ',ILOG,IOB, ** $ 'RIGIDITE',IOB,XVA,' ',ILOG,ri5 ) ** endif if (lagdua.ne.0) > CALL ECCTAB(MTABLE,'ENTIER ',13,XVA,' ',ILOG,IOB, $ 'MAILLAGE',IOB,XVA,' ',ILOG,lagdua) ISUPLO=MTABLE SEGDES RI1,RI2,MTABLE 27 CONTINUE MTABLE=ISUPLO SEGACT MTABLE IF(INSYM.EQ.1) THEN ILUG=.TRUE. ELSE ILUG=.FALSE. ENDIF CALL ECCTAB(MTABLE,'MOT ',4,XVA,'NSYM',ILOG,IOB, $ 'LOGIQUE ',IOB,XVA,' ',ILUG,IOB) if(idepe.ne.0) then * on passe les ideme* a mrem sous forme de listchpo n1=if segini mlchpo,mlchp1 do i=1,if mlchpo.ichpoi(i)=ideme0(1,i) mlchp1.ichpoi(i)=ideme1(1,i) enddo CALL ECCTAB(MTABLE,'ENTIER ',10,XVA,' ',ILOG,IOB, $ 'LISTCHPO',IOB,XVA,' ',ILOG,mlchpo) CALL ECCTAB(MTABLE,'ENTIER ',11,XVA,' ',ILOG,IOB, $ 'LISTCHPO',IOB,XVA,' ',ILOG,mlchp1) * pour mrem on met la derniere raideur condensee. Elle contient les pointeurs pour remonter CALL ECCTAB(MTABLE,'ENTIER ',50,XVA,' ',ILOG,IOB, $ 'RIGIDITE',IOB,XVA,' ',ILOG,ipoiri) endif SEGDES MRIGID DO 26 I=IDEMEM(/1),1,-1 ISECO=IDEMEM(I) CALL ACTOBJ ('CHPOINT ',ISECO,1) CALL ECROBJ ('CHPOINT ',ISECO) 26 CONTINUE SEGSUP IDEMEM CALL ECROBJ ('TABLE ',ISUPLO) SEGINI MTEXTE LTT=8 MTEXT(1:LTT) ='UNILATER' NCART=8 SEGDES MTEXTE CALL ECROBJ('TEXTE',MTEXTE) mrigid=iprigo segdes mrigid norinc=norins RETURN C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ... 30 CONTINUE * il se peut que le dernier chp soit du frottement * on l'enleve car il ne sert a rien si on n'appele pas unilater if (idemem(/1).gt.1.and.idepe.ne.0) then mchpoi=ideme0(idemem(/1),if) segact MCHPOI if (mtypoi.eq.'LX ') idemem(/1)=idemem(/1)-1 endif * frottement SEGDES IDEMEM * write(6,*) ' ichski, igradj,insym ',ichski, igradj,insym * write (6,*) ' imtvid ',imtvid if (imtvid.eq.1) then * matrice vide *** write(6,*) ' attention matrice vide. Système surcontraint ' call erreur(-364) * nsoupo=0 nat=0 segact idemem*mod do i=1,idemem(/1) segini mchpoi mchpoi.ifopoi = ifochs idemem(i)=mchpoi enddo if (noen.eq.0) then call ecrent(0) nbelem=0 nbnn=1 nbsous=0 nbref=0 segini ipt8 ipt8.itypel=1 call actobj('MAILLAGE',ipt8,1) call ecrobj('MAILLAGE',ipt8) endif else * write(6,*) ' appel resou1 -- idemem(1)' * segact idemem * idesec= idemem(1) * call ecchpo(idesec,0) * write(6,*) ' appel resou1 -- ipoiri' * call prrigi ( ipoiri,1) * write(6,*) ' ichski insym ', ichski, insym IF(ICHSKI.EQ.1) CALL RESOU1(IPOIRI,IDEMEM,NOID,NOEN,PREC, > ISTAB,ISOUCI) IF(IGRADJ.EQ.1) CALL GRACO0(IPOIRI,IDEMEM,NOID,NOEN) IF(INSYM .EQ.1) CALL LDMT (IPOIRI,IDEMEM,NOID,NOEN,PREC,ISOUCI) IF(IERR.NE.0) GO TO 5001 endif C C------------------------------------------------------- C LA SOLUTION EST CALCULEE --> ON LA MET EN FORME if (noen.eq.0) then call lirobj('MAILLAGE',ipt8,1,iretou) if (ierr.ne.0) then norinc=norins return endif segact ipt8 call lirent(nben,1,iretou) endif SEGACT IDEMEM*mod N=IDEMEM(/1) do i=1,n mchpoi=idemem(i) * les champs de points qui sortent sont de nature diffuse SEGACT MCHPOI*MOD ifopoi = ifochs NAT = MAX(1,JATTRI(/1)) NSOUPO=IPCHP(/1) SEGADJ MCHPOI JATTRI(1)=1 enddo do 2010 ifois=1,30 segact mrigid mrigid=jrsup if (mrigid.eq.0) goto 2011 if(ierr.ne.0) then norinc=norins return endif call resour(idemem,ideme0,ideme1,mrigid,if,noen,ipt8, > isouci,iverif) if (ierr.ne.0) then norinc=norins return endif if=if-1 2010 continue 2011 continue * * on n'appelle plus verlx car je ne vois pas pourquoi on voudrait que les multiplicateurs de lagrange non éliminés soient nuls * **** call verlx(ri2,iret,mchpo1,noen,ipt8) if (noen.eq.0) then call actobj('MAILLAGE',ipt8,1) call ecrobj('MAILLAGE',ipt8) nben=ipt8.num(/2) call ecrent(nben) endif * do 3 i=1,n iret=idemem(n+1-i) * cas table de liaisons statiques if (itbas.ne.0) then il = n + 1 - i ilo = idnote(il) CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0, & 'TABLE',I1,X1,CHARRE,.true.,ITMOD) if (ierr.ne.0) then norinc=norins return endif c CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN) c CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0, c & 'POINT',0,0.0D0,' ',.TRUE.,IPOIN) CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE', & .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,IRET) else if (ipshpo.gt.0) then mlchpo = ipshpo ichpoi(N+1-I) = iret else CALL ACTOBJ ('CHPOINT ',IRET,1) CALL ECROBJ ('CHPOINT ',IRET) endif 3 CONTINUE c-----fin de boucle sur les solutions C------------------------------------------------------- c MENAGE AVANT DE QUITTER 5001 CONTINUE if (itbas.ne.0) then segdes mtab1 segsup idnote CALL ECROBJ ('TABLE ',itbas) endif if (ipshpo.gt.0) then mlchpo = ipshpo CALL ACTOBJ ('LISTCHPO ',ipshpo,1) CALL ECROBJ ('LISTCHPO ',ipshpo) endif SEGSUP IDEMEM C 5000 CONTINUE norinc=norins END