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    PV090527  23/03/15    21:15:02     11630          
  2.  
  3.       SUBROUTINE RESOU
  4.  
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7. C
  8. C **** CET OPERATEUR SERT A RESOUDRE UN SYSTEME D EQUATIONS LINEAIRES
  9. C ****   CHPOINT = RESOU RIGIDITE CHPOINT
  10. C
  11. C
  12. -INC PPARAM
  13. -INC CCOPTIO
  14.  
  15. -INC SMRIGID
  16. -INC SMTEXTE
  17. -INC SMTABLE
  18. -INC SMCHPOI
  19. -INC SMELEME
  20. -INC SMLCHPO
  21. -INC CCREEL
  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(9)
  29.       CHARACTER*(8) CHARRE1
  30.       CHARACTER*72 CHARRE
  31.       REAL*8 XVA
  32.       LOGICAL ILOG,ILUG,casfimp
  33.       DATA LISM/'NOID','NOUN','ENSE','GRAD','CHOL','STAB','ELIM',
  34.      >'NOST','SOUC'/
  35.       DATA ILOG/.FALSE./
  36. C
  37. C-------------------------------------------------------
  38. c     LECTURE ET INITIALISATION
  39.  
  40. c     LECTURE DES OPTIONS
  41.       XVA=REAL(0.D0)
  42.       IOB=0
  43.  
  44. *  sauvegarde norinc
  45. *     write(6,*) 'norinc norind ',norinc,norind
  46.       norins=norinc
  47.       iverif=1
  48.       ipt8=0
  49.       iunil=0
  50. *  le defaut est de faire une passe d'elimination
  51.       nelim=30
  52. *  experimentalement 2 passes est mieux
  53.       nelim=2
  54.       IMTVID=0
  55.       NOUNIL=0
  56.       NOID=0
  57.       NOEN=1
  58.       IGRADJ = 0
  59.       ICHSKI = 0
  60.       INSYM  = 0
  61.       KIKI=0
  62.       KSYMET = 0
  63.       IPSHPO = 0
  64.       ISTAB=0
  65.       ISOUCI = 0
  66.       ifochs = -99
  67.   5   CONTINUE
  68.       CALL LIRMOT(LISM,9,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) IGRADJ = 1
  73. *      IF (KIKI.EQ.5) ICHSKI = 1
  74. *      IF (KIKI.EQ.4.OR.KIKI.EQ.5) KSYMET = 1
  75.       IF (KIKI.eQ.6) ISTAB=1
  76.       IF (KIKI.eQ.7) then
  77.        call lirent(nelim,1,iretou)
  78.        nelim=min(30,max(0,nelim))
  79.       endif
  80.       IF (KIKI.eQ.8) ISTAB=0
  81.       IF (KIKI.eQ.9) ISOUCI=1
  82.       IF (KIKI.NE.0) GOTO 5
  83.       if (noid.eq.1) iverif=0
  84.       IF(NUCROU.EQ.0) THEN
  85.           ICHSKI=1
  86.       ELSEIF(NUCROU.EQ.1) THEN
  87.           IGRADJ=1
  88.           KSYMET=1
  89.       ENDIF
  90. *      WRITE(6,*) ' nucrou', nucrou
  91. *      IF ( IGRADJ + ICHSKI .EQ. 0 ) ICHSKI = 1
  92.  
  93. c     LECTURE DE LA RIGIDITE
  94.       CALL LIROBJ('RIGIDITE',IPOIRI,1,IRETOU)
  95.       IF(IERR.NE.0) GO TO 5000
  96.       IPRIGO=IPOIRI
  97. C
  98. c     LECTURE DE LA PRECISION
  99.       PREC=REAL(xspeti)
  100.       CALL LIRREE(PREC,0,IRETOU)
  101.       IF(IERR.NE.0) GO TO 5000
  102.  
  103. C     REMPLISSAGE DU 2ND MEMBRE IDEMEM(**) A PARTIR DE ...
  104. c     ... 'CHPOINT'
  105.       SEGINI IDEMEM
  106.    1  CONTINUE
  107.       CALL LIROBJ('CHPOINT ',ISECO,0,IRETOU)
  108.       IF(IRETOU.NE.0) THEN
  109.          IDEMEM(**)=ISECO
  110.          mchpoi=iseco
  111.          segact mchpoi
  112. *         write(6,*) ' extension idemem 1 ',idemem(/1)
  113.          GO TO 1
  114.       ENDIF
  115.  
  116. c     ... 'TABLE DE SOUS-TYPE LIAISONS_STATIQUES'
  117.       CALL LIRTAB('LIAISONS_STATIQUES',ITBAS,0,IRET)
  118.  
  119. c     ... 'LISTCHPO'
  120.       CALL LIROBJ('LISTCHPO',ISECO,0,IRETOU)
  121.       IF(IRETOU.NE.0) THEN
  122.          mlchpo=ISECO
  123.          segact mlchpo
  124.          n1 = ichpoi(/1)
  125.          do iu = 1 , n1
  126.            idemem(**) = ichpoi(iu)
  127. *         write(6,*) ' extension idemem 2 ',idemem(/1)
  128.          enddo
  129.          segdes mlchpo
  130.          segini mlchpo
  131.          ipshpo = mlchpo
  132.       ENDIF
  133.       IF (IERR.NE.0) RETURN
  134.  
  135.       IF (ITBAS.NE.0 .AND. IIMPI.EQ.333) THEN
  136.       WRITE(IOIMP,*) 'on a lu la table des conditions aux limites'
  137.       ENDIF
  138.       if (itbas.ne.0) then
  139.          mtab1 = itbas
  140.          segact mtab1
  141.          ima = mtab1.mlotab - 1
  142.          segini idnote
  143.          im = 0
  144.          segdes mtab1
  145.       else
  146.          goto 90
  147.       endif
  148. *     boucle en cas de résolutions successives avec table
  149.  80   continue
  150.          im = im + 1
  151.          itmod = 0
  152.          ichp0 = 0
  153.          if (im.gt.ima) then
  154.            if (idemem(/1).gt.0) goto 90
  155. *          pas de champs de force
  156.            call erreur(1)
  157.            return
  158.          endif
  159.          CALL ACCTAB(ITBAS,'ENTIER',IM,0.d0,' ',.true.,IP0,
  160.      &                   'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  161.          if (ierr.ne.0) return
  162. c        table itmod trouvee --> on recupere la force
  163.          if (itmod.gt.0) then
  164.            CALL ACCTAB(ITMOD,'MOT',0,0.d0,'FORCE',.true.,IP0,
  165.      &                     'CHPOINT',I1,X1,CHARRE,.true.,ICHP0)
  166.            if (ierr.ne.0) return
  167.            if (ichp0.gt.0) then
  168.              idemem(**) = ichp0
  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.  
  183. C-------------------------------------------------------
  184. c     DEBUT DU TRAVAIL
  185.  
  186.  90   continue
  187.       segini ideme0,ideme1
  188. *  verification pas de blocage en double
  189. ***   call verlag(ipoiri)
  190.       if (ierr.ne.0) return
  191. *        y a t il des matrices de relations non unilaterales
  192.       ipoir0 = ipoiri
  193.       mrigid=ipoiri
  194. C      call prrigi(ipoiri,1)
  195.       segact mrigid
  196.       ifochs = iforig
  197.       nrige= irigel(/1)
  198.       idepe=0
  199. *      write(ioimp,*) 'dans resou mrigid iforig ',mrigid,iforig
  200.  
  201.       nbr = irigel(/2)
  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(jrcond.ne.0) nelim=30
  210.        if (irigel(7,irig).ne.0) then
  211.          insym=1
  212.          ichski=0
  213.        endif
  214. 1000  continue
  215. *  elimination recursive des conditions aux limites
  216. *  on la fait en gradient conjugue ou en appel de unilater
  217.       nfois=nelim-1
  218.        if (igradj.eq.1.or.(iunil.eq.1.and.nounil.eq.0)) nfois=29
  219.       lagdua=0
  220.       imult=1
  221.       icond=idepe
  222.       icondi=(icond*10)/9+1
  223.       if=0
  224.       do ifois=1,nfois
  225.       if(imult.ne.0.and.icond.ne.0.and.(icond*10)/9.lt.icondi.and.
  226.      >      (icondi-icond.gt.0.or.igradj.eq.1)) then
  227.        icondi=icond
  228.        lagdua=-1
  229.        if=if+1
  230.        if(ierr.ne.0) return
  231.        call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
  232.      >       nounil,lagdua,icond,imult,if,imtvid,nelim)
  233. **      write(ioimp,*) ' passe ',if,' condition ',icond,ifois
  234.        if(ierr.ne.0) return
  235.        mrigid=mrigic
  236.       endif
  237.       enddo
  238. *  Si on n'a pas reussi a tout eliminer, on fait encore une passe pour creer lagdua
  239.       lagdua=0
  240.       if (iunil.eq.0.or.nounil.eq.1) then
  241.       if (icond.ne.0) then
  242.        if=if+1
  243.        if(ierr.ne.0) return
  244.         call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
  245.      >        nounil,lagdua,icond,imult,if,imtvid,nelim)
  246. *        write(ioimp,*) ' passe ','finale',' condition ',icond
  247.         if(ierr.ne.0) return
  248.         mrigid=mrigic
  249.       endif
  250.       endif
  251. **     write (ioimp,*) 'nombre de passes',if
  252.       if (idepe.ne.0) noid = 1
  253.       ipoiri=mrigid
  254. *     call prrigi(ipoiri,1)
  255. C-------------------------------------------------------
  256.  
  257. *
  258. * Si au moins une des matrices n'est pas symétrique, on passera
  259. * par le solveur non-symétrique LDMT.
  260. *
  261.       SEGACT MRIGID*MOD
  262.       NRG = IRIGEL(/1)
  263.       NBR = IRIGEL(/2)
  264. C ... Ceci peut arriver si par exemple on extrait la partie
  265. C     symétrique d'une matrice purement antisymétrique ...
  266. *     IF(NBR.EQ.0) THEN
  267. *        SEGDES MRIGID
  268. *        CALL ERREUR(727)
  269. *        RETURN
  270. *     ENDIF
  271. C ... Mais avant on va tester si la normalisation des variables
  272. C     primales et duales a été demandée - ceci entraîne la perte
  273. C     de la symétrie ...
  274.       IF(NORINC.GT.0  .AND. NORIND.GT.0) THEN
  275.          IF(KSYMET.EQ.1) THEN
  276.             CALL ERREUR(19)
  277.             SEGDES,MRIGID
  278.             RETURN
  279.          ENDIF
  280.          INSYM = 1
  281.          IGRADJ = 0
  282.          ICHSKI = 0
  283.          GOTO 15
  284.       ENDIF
  285.  
  286.       IF (NRG.GE.7) THEN
  287. C    ... On teste si la matrice contient des matrices non symétriques ...
  288. C    ... Si OUI, ce n'est pas la peine de faire les autres tests ...
  289.          DO  9 IN = 1,NBR
  290.             IANTI=IRIGEL(7,IN)
  291.             IF(IANTI.GT.0) THEN
  292. C          ... On vérifie si l'utilisateur n'a pas demandé explicitement
  293. C              la résolution par Choleski ou gradient conjugué,
  294. C              si OUI on râle puis on s'en va !!! ...
  295.                IF(KSYMET.EQ.1) THEN
  296.                   CALL ERREUR(1126)
  297.                   SEGDES,MRIGID
  298.                   RETURN
  299.                ENDIF
  300.                IF(NORINC.NE.0.AND.NORIND.EQ.0) THEN
  301. *  on supprime la normalisation au lieu de faire une erreur
  302.                   norinc=0
  303. **                CALL ERREUR(760)
  304. **                SEGDES,MRIGID
  305. **                RETURN
  306.                ENDIF
  307.                INSYM = 1
  308.                IGRADJ = 0
  309.                ICHSKI = 0
  310.                GOTO 15
  311.             ENDIF
  312.   9      CONTINUE
  313.  
  314.  
  315.       ENDIF
  316.  
  317.  15   CONTINUE
  318. C
  319. C  ****  ON S'ASSURE QU'IL N'Y A PAS D'APPUIS UNILATERAUX
  320. C
  321.       if (iunil.eq.0) goto 30
  322.       IF(IRIGEL(/1).LT.6) GO TO 30
  323.       IF (NOUNIL.EQ.1) GOTO 30
  324.    21 CONTINUE
  325. C
  326. C  **** EXISTENCE DES APPUIS UNILATERAUX
  327. C  **** SI ON EST DEJA PASSE DANS LES APPUIS UNILATERAUX
  328. C       ISUPEQ POINTE SUR UNE TABLE CONTENANT LES DONNEES A PASSER
  329. C       A LA PROCEDURE UNILATER
  330. C
  331.       ISUPLO=ISUPEQ
  332.       IF (ISUPLO.NE.0) GOTO 27
  333.       NNOR=0
  334.       DO 22 I=1,IRIGEL(/2)
  335.       IF(IRIGEL(6,I).EQ.0) NNOR=NNOR+1
  336.   22  CONTINUE
  337.       IF(NNOR.EQ.0)  THEN
  338.                        CALL ERREUR(312)
  339.                        norinc=norins
  340.                        RETURN
  341.                      ENDIF
  342.       NRIGE=IRIGEL(/1)
  343.       NRIGEL=NNOR
  344.       SEGINI RI1
  345.       RI1.IFORIG = IFORIG
  346. c*      RI1.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
  347.       RI1.MTYMAT = '        '
  348.       NRIGEL=IRIGEL(/2)-NNOR
  349.       SEGINI RI2
  350.       RI2.IFORIG = IFORIG
  351. c*      RI2.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
  352.       RI2.MTYMAT = '        '
  353.       II1=0
  354.       II2=0
  355.       DO 23 I=1,IRIGEL(/2)
  356.       IF(IRIGEL(6,I).NE.0) THEN
  357.                                 RI3=RI2
  358.                                 II2=II2+1
  359.                                 II=II2
  360.                            ELSE
  361.                                 RI3=RI1
  362.                                 II1=II1+1
  363.                                 II=II1
  364.                            ENDIF
  365.       DO 24 J=1,NRIGE
  366.       RI3.IRIGEL(J,II) = IRIGEL(J,I)
  367.    24 CONTINUE
  368.       RI3.COERIG(II)=COERIG(I)
  369.    23 CONTINUE
  370. *  RI1 raideur sans condition unilaterale
  371. *  RI2 conditions unilaterales
  372.       CALL CRTABL(MTABLE)
  373.       ISUPEQ=MTABLE
  374. * il faut aussi mettre isupeq dans la raideur initiale
  375.       if (jrsup.ne.0) mrigid=jrsup
  376.       segact mrigid
  377.       iri1s=jrelim
  378.       iri2s=mrigid
  379.       MRIGID=IPRIGO
  380.       SEGACT MRIGID*MOD
  381.       ISUPEQ=MTABLE
  382.       if (idepe.ne.0) then
  383. *  il faut extraire de la matrice initiale (ipoir0) les conditions unilaterales
  384.       mrigid=iri2s
  385.       segact mrigid
  386.       nrigel=0
  387.       do 40 i=1,irigel(/2)
  388.        if (irigel(6,i).eq.0) nrigel=nrigel+1
  389.   40  continue
  390.       if (nrigel.eq.0) call erreur(312)
  391.       if (ierr.ne.0) then
  392.             norinc=norins
  393.             return
  394.       endif
  395.       nrige=irigel(/1)
  396.       segini ri4
  397.       ri4.iforig = iforig
  398.       ri4.mtymat = mtymat
  399.       ii1=0
  400.       nrigel=irigel(/2)-nrigel
  401.       segini ri5
  402.       ri5.iforig = iforig
  403.       ri5.mtymat = mtymat
  404.       ii2=0
  405.       do 41 i=1,irigel(/2)
  406.        if (irigel(6,i).ne.0) goto 42
  407.        ii1=ii1+1
  408.        do j=1,nrige
  409.         ri4.irigel(j,ii1)=irigel(j,i)
  410.        enddo
  411.        ri4.coerig(ii1)=coerig(i)
  412.        goto 41
  413.   42  continue
  414.        ii2=ii2+1
  415.        do j=1,nrige
  416.         ri5.irigel(j,ii2)=irigel(j,i)
  417.        enddo
  418.        ri5.coerig(ii2)=coerig(i)
  419.   41  continue
  420.       segdes mrigid,ri4
  421.       endif
  422.       ri3=iri1s
  423. *     segact ri1,ri2,ri3,ri4,ri5
  424.       CALL ECCTAB(MTABLE,'ENTIER  ',1,XVA,' ',ILOG,IOB,
  425.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI1)
  426.       CALL ECCTAB(MTABLE,'ENTIER  ',2,XVA,' ',ILOG,IOB,
  427.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI2)
  428.       CALL ECCTAB(MTABLE,'ENTIER  ',3,XVA,' ',ILOG,IOB,
  429.      $                   'LOGIQUE ',IOB,XVA,' ',ILOG,IOB)
  430. **      if(idepe.ne.0) then
  431. **    CALL ECCTAB(MTABLE,'ENTIER  ',8,XVA,' ',ILOG,IOB,
  432. **   $                   'RIGIDITE',IOB,XVA,' ',ILOG,iri1s)
  433. **    CALL ECCTAB(MTABLE,'ENTIER  ',9,XVA,' ',ILOG,IOB,
  434. **   $                   'RIGIDITE',IOB,XVA,' ',ILOG,ri4  )
  435. **    CALL ECCTAB(MTABLE,'ENTIER  ',12,XVA,' ',ILOG,IOB,
  436. **   $                   'RIGIDITE',IOB,XVA,' ',ILOG,ri5  )
  437. **      endif
  438.         if (lagdua.ne.0)
  439.      > CALL ECCTAB(MTABLE,'ENTIER  ',13,XVA,' ',ILOG,IOB,
  440.      $                   'MAILLAGE',IOB,XVA,' ',ILOG,lagdua)
  441.       ISUPLO=MTABLE
  442.       SEGDES RI1,RI2,MTABLE
  443.   27  CONTINUE
  444.       MTABLE=ISUPLO
  445.       SEGACT MTABLE
  446.       IF(INSYM.EQ.1) THEN
  447.          ILUG=.TRUE.
  448.       ELSE
  449.          ILUG=.FALSE.
  450.       ENDIF
  451.       CALL ECCTAB(MTABLE,'MOT     ',4,XVA,'NSYM',ILOG,IOB,
  452.      $                   'LOGIQUE ',IOB,XVA,' ',ILUG,IOB)
  453.       if(idepe.ne.0) then
  454. *  on passe les ideme* a mrem sous forme de listchpo
  455.         n1=if
  456.         segini mlchpo,mlchp1
  457.         do i=1,if
  458.           mlchpo.ichpoi(i)=ideme0(1,i)
  459.           mlchp1.ichpoi(i)=ideme1(1,i)
  460.         enddo
  461.         CALL ECCTAB(MTABLE,'ENTIER  ',10,XVA,' ',ILOG,IOB,
  462.      $                   'LISTCHPO',IOB,XVA,' ',ILOG,mlchpo)
  463.         CALL ECCTAB(MTABLE,'ENTIER  ',11,XVA,' ',ILOG,IOB,
  464.      $                   'LISTCHPO',IOB,XVA,' ',ILOG,mlchp1)
  465. *  pour mrem on met la derniere raideur condensee. Elle contient les pointeurs pour remonter
  466.         CALL ECCTAB(MTABLE,'ENTIER  ',50,XVA,' ',ILOG,IOB,
  467.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,ipoiri)
  468.       endif
  469.       SEGDES MRIGID
  470.       DO 26 I=IDEMEM(/1),1,-1
  471.          ISECO=IDEMEM(I)
  472.          CALL ACTOBJ ('CHPOINT ',ISECO,1)
  473.          CALL ECROBJ ('CHPOINT ',ISECO)
  474.   26  CONTINUE
  475.       SEGSUP IDEMEM
  476.       CALL ECROBJ ('TABLE   ',ISUPLO)
  477.       SEGINI MTEXTE
  478.       LTT=8
  479.       MTEXT(1:LTT) ='UNILATER'
  480.       NCART=8
  481.       SEGDES MTEXTE
  482.       CALL ECROBJ('TEXTE',MTEXTE)
  483.       mrigid=iprigo
  484.       segdes mrigid
  485.             norinc=norins
  486.       RETURN
  487.  
  488. C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ...
  489.    30 CONTINUE
  490. * il se peut que le dernier chp soit du frottement
  491. * on l'enleve car il ne sert a rien si on n'appele pas unilater
  492.         if (idemem(/1).gt.1.and.idepe.ne.0) then
  493.          mchpoi=ideme0(idemem(/1),if)
  494.          segact MCHPOI
  495.          if (mtypoi.eq.'LX  ') idemem(/1)=idemem(/1)-1
  496.         endif
  497. *  frottement
  498.       SEGDES IDEMEM
  499. *      write(6,*) ' ichski, igradj,insym ',ichski, igradj,insym
  500. *     write (6,*) ' imtvid ',imtvid
  501.       if (imtvid.eq.1) then
  502. *  matrice vide
  503. ***     write(6,*) ' attention matrice vide. Système surcontraint '
  504.         call erreur(-364)
  505. *
  506.         nsoupo=0
  507.         nat=0
  508.         segact idemem*mod
  509.  
  510.         do i=1,idemem(/1)
  511.          segini mchpoi
  512.          mchpoi.ifopoi = ifochs
  513.          idemem(i)=mchpoi
  514.         enddo
  515.         if (noen.eq.0) then
  516.           call ecrent(0)
  517.           nbelem=0
  518.           nbnn=1
  519.           nbsous=0
  520.           nbref=0
  521.           segini ipt8
  522.           ipt8.itypel=1
  523.           call actobj('MAILLAGE',ipt8,1)
  524.           call ecrobj('MAILLAGE',ipt8)
  525.         endif
  526.       else
  527. *       write(6,*) ' appel resou1 --  idemem(1)'
  528. *       segact idemem
  529. *       idesec= idemem(1)
  530. *       call ecchpo(idesec,0)
  531. *       write(6,*) ' appel resou1 -- ipoiri'
  532. *       call prrigi ( ipoiri,1)
  533. *       write(6,*) ' ichski  insym  ', ichski, insym
  534.        IF(ICHSKI.EQ.1) CALL RESOU1(IPOIRI,IDEMEM,NOID,NOEN,PREC,
  535.      >        ISTAB,ISOUCI)
  536.        IF(IGRADJ.EQ.1) CALL GRACO0(IPOIRI,IDEMEM,NOID,NOEN)
  537.        IF(INSYM .EQ.1) CALL LDMT  (IPOIRI,IDEMEM,NOID,NOEN,PREC,ISOUCI)
  538.        IF(IERR.NE.0) GO TO 5001
  539.       endif
  540. C
  541. C-------------------------------------------------------
  542. C     LA SOLUTION EST CALCULEE --> ON LA MET EN FORME
  543.  
  544.        if (noen.eq.0) then
  545.          call lirobj('MAILLAGE',ipt8,1,iretou)
  546.          if (ierr.ne.0) then
  547.             norinc=norins
  548.             return
  549.          endif
  550.          segact ipt8
  551.          call lirent(nben,1,iretou)
  552.        endif
  553.  
  554.       SEGACT IDEMEM*mod
  555.       N=IDEMEM(/1)
  556.       do i=1,n
  557.         mchpoi=idemem(i)
  558. *     les champs de points qui sortent sont de nature diffuse
  559.          SEGACT MCHPOI*MOD
  560.          ifopoi = ifochs
  561.          NAT = MAX(1,JATTRI(/1))
  562.          NSOUPO=IPCHP(/1)
  563.          SEGADJ MCHPOI
  564.          JATTRI(1)=1
  565.       enddo
  566.  
  567.       do 2010 ifois=1,30
  568.       segact mrigid
  569.       mrigid=jrsup
  570.       if (mrigid.eq.0) goto 2011
  571.       if(ierr.ne.0) then
  572.             norinc=norins
  573.             return
  574.          endif
  575.       call resour(idemem,ideme0,ideme1,mrigid,if,noen,ipt8,
  576.      >       isouci,iverif)
  577.       if (ierr.ne.0) then
  578.             norinc=norins
  579.             return
  580.          endif
  581.       if=if-1
  582.  2010 continue
  583.  2011 continue
  584. *
  585. *  on n'appelle plus verlx car je ne vois pas pourquoi on voudrait que les multiplicateurs de lagrange non éliminés soient nuls
  586. *
  587. ****   call verlx(ri2,iret,mchpo1,noen,ipt8)
  588.        if (noen.eq.0) then
  589.          call actobj('MAILLAGE',ipt8,1)
  590.          call ecrobj('MAILLAGE',ipt8)
  591.          nben=ipt8.num(/2)
  592.          call ecrent(nben)
  593.        endif
  594. *
  595.       do 3 i=1,n
  596.        iret=idemem(n+1-i)
  597.  
  598. *     cas table de liaisons statiques
  599.       if (itbas.ne.0) then
  600.            il = n + 1 - i
  601.            ilo = idnote(il)
  602.          CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0,
  603.      &                   'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  604.            if (ierr.ne.0) then
  605.             norinc=norins
  606.             return
  607.          endif
  608. c          CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
  609. c          CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  610. c      &        'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  611.  
  612.          CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE',
  613.      &        .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,IRET)
  614.  
  615.       else if (ipshpo.gt.0)  then
  616.          mlchpo = ipshpo
  617.          ichpoi(N+1-I) = iret
  618.       else
  619.          CALL ACTOBJ ('CHPOINT ',IRET,1)
  620.          CALL ECROBJ ('CHPOINT ',IRET)
  621.       endif
  622.  
  623.    3  CONTINUE
  624. c-----fin de boucle sur les solutions
  625.  
  626.  
  627. C-------------------------------------------------------
  628. c     MENAGE AVANT DE QUITTER
  629.  
  630.  5001 CONTINUE
  631.       if (itbas.ne.0) then
  632.          segdes mtab1
  633.          segsup idnote
  634.          CALL ECROBJ ('TABLE   ',itbas)
  635.       endif
  636.  
  637.       if (ipshpo.gt.0) then
  638.         mlchpo = ipshpo
  639.         CALL ACTOBJ ('LISTCHPO ',ipshpo,1)
  640.         CALL ECROBJ ('LISTCHPO ',ipshpo)
  641.       endif
  642.       SEGSUP IDEMEM
  643. C
  644.  5000 CONTINUE
  645.       norinc=norins
  646.       END
  647.  
  648.  
  649.  
  650.  
  651.  
  652.  
  653.  

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