Télécharger la plaquette Cast3M


fr en

Télécharger resou.eso

Retour à la liste

Numérotation des lignes : 

resou

  1. C RESOU     SOURCE    MB234859  26/06/25    21:15:21     12580          
  2.  
  3.       SUBROUTINE RESOU
  4. C----------------------------------------------------------------------
  5. C **** CET OPERATEUR SERT A RESOUDRE UN SYSTEME D EQUATIONS LINEAIRES
  6. C ****   CHPOINT = RESOU RIGIDITE CHPOINT
  7. C----------------------------------------------------------------------
  8.       IMPLICIT INTEGER(I-N)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10. C
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC SMRIGID
  14. -INC SMCOORD
  15. -INC SMTEXTE
  16. -INC SMTABLE
  17. -INC SMCHPOI
  18. -INC SMELEME
  19. -INC SMLCHPO
  20. -INC CCREEL
  21. C
  22.       PARAMETER(ZERO=0.D0)
  23.       SEGMENT IDEMEM(0)
  24.       segment ideme0(idemem(/1),30)
  25.       segment ideme1(idemem(/1),30)
  26.       segment idnote(0)
  27. C
  28.       CHARACTER*4  LISM(7)
  29.       CHARACTER*72 CHARRE
  30.       REAL*8 XVA
  31.       LOGICAL ILOG,ILUG,casfimp,notok,bdblx
  32.       DATA LISM/'NOID','NOUN','ENSE','STAB','ELIM','NOST','SOUC'/
  33.       DATA ILOG/.FALSE./
  34. C
  35.       XVA=REAL(0.D0)
  36.       IOB=0
  37. C
  38. C     Sauvegarde norinc
  39. *     write(6,*) 'norinc norind ',norinc,norind
  40.       norins=norinc
  41.       ipt8=0
  42.       iunil=0
  43.       IMTVID=0
  44.       IGRADJ=0
  45.       IF (NUCROU.EQ.1) IGRADJ=1
  46.       INSYM=0
  47. C     Tester si la normalisation des variables primales et duales a
  48. C     été demandée - ceci entraîne la perte de la symétrie ...
  49.       IF (NORINC.GT.0  .AND. NORIND.GT.0) THEN
  50.         IF (IGRADJ.EQ.1) THEN
  51.           CALL ERREUR(19)
  52.           SEGDES,MRIGID
  53.           RETURN
  54.         ENDIF
  55.         INSYM=1
  56.       ENDIF
  57.       IPSHPO = 0
  58.       ifochs = -99
  59.       NDEMEM=0
  60. C----------------------------------------------------------------------
  61. C     LECTURE DES ARGUMENTS DE L'OPERATEUR RESOU
  62. C----------------------------------------------------------------------
  63. C     LECTURE DES OPTIONS/MOTS-CLES
  64.       KIKI=0
  65.       NOID=0
  66.       NOUNIL=0
  67.       NOEN=1
  68.       ISTAB=0
  69.       NELIM=30
  70.       ISOUCI=0
  71.   5   CONTINUE
  72.       CALL LIRMOT(LISM,7,KIKI,0)
  73.       IF (KIKI.EQ.1) NOID=1
  74.       IF (KIKI.EQ.2) NOUNIL=1
  75.       IF (KIKI.EQ.3) NOEN=0
  76.       IF (KIKI.EQ.4) ISTAB=1
  77.       IF (KIKI.EQ.5) THEN
  78.         CALL LIRENT(NELIM,1,IRETOU)
  79.         NELIM=MIN(30,MAX(0,NELIM))
  80.       ENDIF
  81.       IF (KIKI.EQ.6) ISTAB=0
  82.       IF (KIKI.EQ.7) ISOUCI=1
  83.       IF (KIKI.NE.0) GOTO 5
  84. C
  85.       iverif=1
  86.       if (noid.eq.1) iverif=0
  87. C
  88. C     LECTURE DE LA RIGIDITE
  89.       CALL LIROBJ('RIGIDITE',IPOIRI,1,IRETOU)
  90.       IF(IERR.NE.0) GO TO 5000
  91.       IPRIGO=IPOIRI
  92. C
  93. C     LECTURE DE LA PRECISION
  94.       PREC=REAL(xspeti)
  95.       CALL LIRREE(PREC,0,IRETOU)
  96.       IF(IERR.NE.0) GO TO 5000
  97. C
  98. C     REMPLISSAGE DU 2ND MEMBRE IDEMEM(**) A PARTIR DE
  99. C     ... CHPOINT
  100.       SEGINI IDEMEM
  101.    1  CONTINUE
  102.       CALL LIROBJ('CHPOINT ',ISECO,0,IRETOU)
  103.       IF(IRETOU.NE.0) THEN
  104.         IDEMEM(**)=ISECO
  105.         NDEMEM=NDEMEM+1
  106.         mchpoi=iseco
  107.         segact mchpoi
  108.         if (mchpoi.jattri(/1).ge.1) then
  109.           notok=(mchpoi.jattri(1).ne.2)
  110.         else
  111.           notok=.true.
  112.         endif
  113.         if (notok) then
  114.           INTERR=mchpoi
  115.           CALL ERREUR(-387)
  116.         endif
  117.         GOTO 1
  118.       ENDIF
  119. C
  120. C     ... LISTCHPO
  121.       CALL LIROBJ('LISTCHPO',ISECO,0,IRETOU)
  122.       IF(IRETOU.NE.0) THEN
  123.         mlchpo=ISECO
  124.         segact mlchpo
  125.         NDEMEM = ichpoi(/1)
  126.         do iu = 1,NDEMEM
  127.           idemem(**) = ichpoi(iu)
  128. *         write(6,*) ' extension idemem 2 ',idemem(/1)
  129.         enddo
  130.         segdes mlchpo
  131.         n1 = ndemem
  132.         segini mlchpo
  133.         ipshpo = mlchpo
  134.       ENDIF
  135.       IF (IERR.NE.0) RETURN
  136. C
  137. C     ... TABLE DE SOUS-TYPE LIAISONS_STATIQUES
  138.       CALL LIRTAB('LIAISONS_STATIQUES',ITBAS,0,IRET)
  139.       IF (ITBAS.EQ.0) GOTO 90
  140. C
  141.       IF (IIMPI.EQ.333) THEN
  142.         WRITE(IOIMP,*) 'on a lu la table des conditions aux limites'
  143.       ENDIF
  144.       mtab1 = itbas
  145.       segact mtab1
  146.       ima = mtab1.mlotab - 1
  147.       segini idnote
  148.       im = 0
  149.       segdes mtab1
  150.  80   CONTINUE
  151.       im = im + 1
  152.       itmod = 0
  153.       ichp0 = 0
  154.       if (im.gt.ima) then
  155.         IF (NDEMEM.GT.0) GOTO 90
  156. *       pas de champs de force
  157.         CALL ERREUR(1)
  158.         return
  159.       endif
  160.       CALL ACCTAB(ITBAS,'ENTIER',IM,0.d0,' ',.true.,IP0,
  161.      &                 'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  162.       if (ierr.ne.0) return
  163. c     table itmod trouvee --> on recupere la force
  164.       if (itmod.gt.0) then
  165.         CALL ACCTAB(ITMOD,'MOT',0,0.d0,'FORCE',.true.,IP0,
  166.      &                  'CHPOINT',I1,X1,CHARRE,.true.,ICHP0)
  167.         if (ierr.ne.0) return
  168.         if (ichp0.gt.0) then
  169.           idemem(**) = ichp0
  170.           NDEMEM=NDEMEM+1
  171. *       write(6,*) ' extension idemem 3 ',idemem(/1)
  172.           idnote(**) = im
  173.         else
  174.           call erreur(1)
  175.           return
  176.         endif
  177. c       on cree le point repere ici
  178.         CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
  179.         CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  180.      &        'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  181.       endif
  182.       GOTO 80
  183.       IF (IERR.NE.0) RETURN
  184. C----------------------------------------------------------------------
  185. C     DEBUT DU TRAVAIL
  186.  90   CONTINUE
  187.       SEGINI IDEME0,IDEME1
  188. C
  189. C     Verifier qu'il n'y a pas de blocage en double
  190. ***   call verlag(ipoiri)
  191.       if (ierr.ne.0) return
  192. *     y a t il des matrices de relations non unilaterales
  193.       mrigid=ipoiri
  194. C     call prrigi(ipoiri,1)
  195.       segact mrigid
  196.       ifochs=iforig
  197.       idepe=0
  198. *     write(ioimp,*) 'dans resou mrigid iforig ',mrigid,iforig
  199. C
  200.       nbr = irigel(/2)
  201.       if (jrcond.ne.0) nelim=30
  202.       do 1000 irig = 1,nbr
  203.         meleme=irigel(1,irig)
  204.         segact meleme
  205.         if ((irigel(6,irig).eq.0.or.nounil.eq.1).and.itypel.eq.22)
  206.      >    idepe=idepe+num(/2)
  207.         if (irigel(6,irig).ne.0) iunil=1
  208.         if (irigel(6,irig).eq.2) nelim=30
  209.         if (irigel(7,irig).ne.0) insym=1
  210. 1000  continue
  211. C
  212. C     Elimination recursive des conditions aux limites
  213. *  on la fait en gradient conjugue ou en appel de unilater
  214.       if (igradj.eq.1.or.(iunil.eq.1.and.nounil.eq.0)) nelim=30
  215.       nfois=nelim-1
  216.       bdblx=.false.
  217.       imult=1
  218.       icond=idepe
  219.       icondi=(icond*10)/9+1
  220.       if=0
  221.       do ifois=1,nfois
  222.         if(imult.ne.0.and.icond.ne.0.and.(icond*10)/9.lt.icondi.and.
  223.      >      (icondi-icond.gt.0.or.igradj.eq.1)) then
  224.           icondi=icond
  225.           if=if+1
  226.           if(ierr.ne.0) return
  227.           call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
  228.      >                nounil,bdblx,icond,imult,if,imtvid,nelim)
  229. C         write(ioimp,*) ' passe ',if,' condition ',icond,ifois
  230.           if(ierr.ne.0) return
  231.           mrigid=mrigic
  232.         endif
  233.       enddo
  234. C
  235. C     S'il reste des conditions : dedoubler les mult de Lagrange restants
  236. C     -> nouvel appel pour creer lagdua et adapter les seconds membres
  237.       if (iunil.eq.0.or.nounil.eq.1) then
  238.         if (icond.ne.0) then
  239.           if=if+1
  240.           bdblx=.true.
  241.           if(ierr.ne.0) return
  242.           call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
  243.      >                nounil,bdblx,icond,imult,if,imtvid,nelim)
  244. *         write(ioimp,*) ' passe ','finale',' condition ',icond
  245.           if(ierr.ne.0) return
  246.           mrigid=mrigic
  247.         endif
  248.       endif
  249. *     write (ioimp,*) 'nombre de passes',if,mrigid.imlag
  250.       if (idepe.ne.0) noid=1
  251. C
  252. C     Rigidite reduite aux inconnues independantes
  253.       ipoiri=mrigid
  254. *     call prrigi(ipoiri,1)
  255. C-------------------------------------------------------
  256. *
  257. * Si au moins une des matrices n'est pas symétrique, on passera
  258. * par le solveur non-symétrique LDMT.
  259. *
  260.       SEGACT MRIGID*MOD
  261. CCC   IF (INSYM.EQ.1) GOTO 15
  262. CCC   NBR = IRIGEL(/2)
  263. C ... Ceci peut arriver si par exemple on extrait la partie
  264. C     symétrique d'une matrice purement antisymétrique ...
  265. *     IF(NBR.EQ.0) THEN
  266. *        SEGDES MRIGID
  267. *        CALL ERREUR(727)
  268. *        RETURN
  269. *     ENDIF
  270. C
  271. C     DO 9 IN = 1,NBR
  272. C       IF (IRIGEL(7,IN).GT.0) THEN
  273. C         INSYM=1
  274. C         GOTO 15
  275. C       ENDIF
  276. C 9   CONTINUE
  277. C15   CONTINUE
  278.  
  279.       IF (INSYM.EQ.1) THEN
  280. C      ... On vérifie si l'utilisateur n'a pas demandé explicitement
  281. C          la résolution par Choleski ou gradient conjugué,
  282. C          si OUI on râle puis on s'en va !!! ...
  283.         IF (IGRADJ.EQ.1) THEN
  284.           CALL ERREUR(1126)
  285.           SEGDES,MRIGID
  286.           RETURN
  287.         ENDIF
  288.         IF (NORINC.NE.0.AND.NORIND.EQ.0) THEN
  289. *  on supprime la normalisation au lieu de faire une erreur
  290.           norinc=0
  291. **        CALL ERREUR(760)
  292. **        SEGDES,MRIGID
  293. **        RETURN
  294.         ENDIF
  295.       ENDIF
  296. C
  297.       IF (IUNIL.EQ.0.OR.NOUNIL.EQ.1) GOTO 30
  298. C
  299. C  **** EXISTENCE DES APPUIS UNILATERAUX
  300. C  **** SI ON EST DEJA PASSE DANS LES APPUIS UNILATERAUX
  301. C       ISUPEQ POINTE SUR UNE TABLE CONTENANT LES DONNEES A PASSER
  302. C       A LA PROCEDURE UNILATER
  303. C
  304.       ISUPLO=ISUPEQ
  305.       IF (ISUPLO.NE.0) GOTO 27
  306. C
  307. C     Distinguer ce qui est unilateral (RI2) du reste (RI1)
  308.       NNOR=0
  309.       DO 22 I=1,IRIGEL(/2)
  310.         IF(IRIGEL(6,I).EQ.0) NNOR=NNOR+1
  311.   22  CONTINUE
  312.       IF(NNOR.EQ.0)  THEN
  313.         CALL ERREUR(312)
  314.         GOTO 5000
  315.       ENDIF
  316. C
  317.       NRIGE=IRIGEL(/1)
  318.       NRIGEL=NNOR
  319.       SEGINI RI1
  320.       RI1.IFORIG = IFORIG
  321. c*      RI1.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
  322.       RI1.MTYMAT = '        '
  323.       NRIGEL=IRIGEL(/2)-NNOR
  324.       SEGINI RI2
  325.       RI2.IFORIG = IFORIG
  326. c*      RI2.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
  327.       RI2.MTYMAT = '        '
  328.       II1=0
  329.       II2=0
  330.       DO 23 I=1,IRIGEL(/2)
  331.         IF(IRIGEL(6,I).NE.0) THEN
  332.           RI3=RI2
  333.           II2=II2+1
  334.           II=II2
  335.         ELSE
  336.           RI3=RI1
  337.           II1=II1+1
  338.           II=II1
  339.         ENDIF
  340.         DO 24 J=1,NRIGE
  341.           RI3.IRIGEL(J,II) = IRIGEL(J,I)
  342.    24   CONTINUE
  343.         RI3.COERIG(II)=COERIG(I)
  344.    23 CONTINUE
  345. C
  346. C     Creation de la table a transmettre a UNILATER (stockee dans la
  347. C     raideur initiale)
  348.       CALL CRTABL(MTABLE)
  349.       ISUPEQ=MTABLE
  350.       MRIGID=IPRIGO
  351.       SEGACT MRIGID*MOD
  352.       ISUPEQ=MTABLE
  353.       CALL ECCTAB(MTABLE,'ENTIER  ',1  ,XVA,' ',ILOG,IOB,
  354.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI1)
  355.       CALL ECCTAB(MTABLE,'ENTIER  ',2  ,XVA,' ',ILOG,IOB,
  356.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI2)
  357.       CALL ECCTAB(MTABLE,'ENTIER  ',3  ,XVA,' ',ILOG,IOB,
  358.      $                   'LOGIQUE ',IOB,XVA,' ',ILOG,IOB)
  359.  
  360.       ISUPLO=MTABLE
  361.       SEGDES RI1,RI2,MTABLE
  362. C
  363. C     Recuperer la table a transmettre a UNILATER
  364.   27  CONTINUE
  365.       MTABLE=ISUPLO
  366.       SEGACT MTABLE
  367.       ILUG=(INSYM.EQ.1)
  368. C
  369.       CALL ECCTAB(MTABLE,'MOT     ',4  ,XVA,'NSYM',ILOG,IOB,
  370.      $                   'LOGIQUE ',IOB,XVA,' '   ,ILUG,IOB)
  371.       if(idepe.ne.0) then
  372. *  on passe les ideme* a mrem sous forme de listchpo
  373.         n1=if
  374.         segini mlchpo,mlchp1
  375.         do i=1,if
  376.           mlchpo.ichpoi(i)=ideme0(1,i)
  377.           mlchp1.ichpoi(i)=ideme1(1,i)
  378.         enddo
  379.         CALL ECCTAB(MTABLE,'ENTIER  ',10 ,XVA,' ',ILOG,IOB,
  380.      $                     'LISTCHPO',IOB,XVA,' ',ILOG,mlchpo)
  381.         CALL ECCTAB(MTABLE,'ENTIER  ',11 ,XVA,' ',ILOG,IOB,
  382.      $                     'LISTCHPO',IOB,XVA,' ',ILOG,mlchp1)
  383. *  pour mrem on met la derniere raideur condensee. Elle contient les pointeurs pour remonter
  384.         CALL ECCTAB(MTABLE,'ENTIER  ',50 ,XVA,' ',ILOG,IOB,
  385.      $                     'RIGIDITE',IOB,XVA,' ',ILOG,ipoiri)
  386.       endif
  387.       SEGDES MRIGID
  388.       DO 26 I=NDEMEM,1,-1
  389.         ISECO=IDEMEM(I)
  390.         CALL ACTOBJ ('CHPOINT ',ISECO,1)
  391.         CALL ECROBJ ('CHPOINT ',ISECO)
  392.   26  CONTINUE
  393.       SEGSUP IDEMEM
  394.       CALL ECROBJ ('TABLE   ',ISUPLO)
  395.       SEGINI MTEXTE
  396.       LTT=8
  397.       MTEXT(1:LTT) ='UNILATER'
  398.       NCART=8
  399.       SEGDES MTEXTE
  400.       CALL ECROBJ('TEXTE',MTEXTE)
  401.       mrigid=iprigo
  402.       segdes mrigid
  403.       GOTO 5000
  404. C
  405. C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ...
  406.    30 CONTINUE
  407. * il se peut que le dernier chp soit du frottement
  408. * on l'enleve car il ne sert a rien si on n'appele pas unilater
  409.       if (NDEMEM.gt.1.and.idepe.ne.0) then
  410.         mchpoi=ideme0(NDEMEM,if)
  411.         segact MCHPOI
  412.         if (mtypoi.eq.'LX  ') NDEMEM=NDEMEM-1
  413.       endif
  414. C
  415. C     Cas matrice vide
  416.       IF (IMTVID.EQ.1) THEN
  417. ***     write(6,*) ' attention matrice vide. Système surcontraint '
  418.         if (nounil.eq.0) call erreur(-364)
  419. *
  420.         nsoupo=0
  421.         nat=0
  422.         segact idemem*mod
  423.  
  424.         do i=1,idemem(/1)
  425.           segini mchpoi
  426.           mchpoi.ifopoi = ifochs
  427.           idemem(i)=mchpoi
  428.         enddo
  429.         if (noen.eq.0) then
  430.           nat=2
  431.           nsoupo=0
  432.           segini mchpo4
  433.           call ecrobj('CHPOINT ',mchpo4)
  434.           call ecrent(0)
  435.         endif
  436.       ELSE
  437.         CALL RESOU1(IPOIRI,IDEMEM,
  438.      &              NOID,NOEN,PREC,ISTAB,ISOUCI,INSYM,IGRADJ)
  439.         IF(IERR.NE.0) GO TO 5000
  440.       ENDIF
  441. C
  442. C     -----------------------------------------------------------------
  443. C     Reintroduire les inconnues eliminees
  444.       do 2010 ifois=1,30
  445.         segact mrigid
  446.         mrigid=jrsup
  447.         if (mrigid.eq.0) goto 2011
  448.         if(ierr.ne.0) GOTO 5000
  449.         call resour(idemem,ideme0,ideme1,mrigid,if,ipt8,isouci,iverif)
  450.         if(ierr.ne.0) GOTO 5000
  451.         if=if-1
  452.  2010 continue
  453.  2011 continue
  454. C
  455. C     -----------------------------------------------------------------
  456. C     ECRITURE DES OBJETS RESULTATS
  457.       SEGACT IDEMEM
  458.       do 3 i=1,NDEMEM
  459.        il = NDEMEM + 1 - i
  460.        iret=idemem(il)
  461. C
  462.        mchpoi=iret
  463.        segact mchpoi*mod
  464.        mchpoi.jattri(1)=1
  465. C
  466.       if (itbas.ne.0) then
  467.          ilo = idnote(il)
  468.          CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0,
  469.      &                   'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  470.          if (ierr.ne.0) GOTO 5000
  471.  
  472.          CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE',
  473.      &        .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,IRET)
  474.  
  475.          if (i.EQ.NDEMEM) then
  476.            segdes mtab1
  477.            segsup idnote
  478.            CALL ECROBJ ('TABLE   ',itbas)
  479.         endif
  480.       else if (ipshpo.gt.0) then
  481.          mlchpo = ipshpo
  482.          ichpoi(il) = iret
  483.          if (i.EQ.NDEMEM) then
  484.           mlchpo = ipshpo
  485.           CALL ACTOBJ ('LISTCHPO ',ipshpo,1)
  486.           CALL ECROBJ ('LISTCHPO ',ipshpo)
  487.         endif
  488.       else
  489.          CALL ACTOBJ ('CHPOINT ',IRET,1)
  490.          CALL ECROBJ ('CHPOINT ',IRET)
  491.       endif
  492.  
  493.    3  CONTINUE
  494.       SEGSUP IDEMEM
  495. C
  496.  5000 CONTINUE
  497.       norinc=norins
  498.       END
  499.  
  500.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales