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

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