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  24/11/11    21:15:04     12068          
  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.  
  442.       ISUPLO=MTABLE
  443.       SEGDES RI1,RI2,MTABLE
  444.   27  CONTINUE
  445.       MTABLE=ISUPLO
  446.       SEGACT MTABLE
  447.       IF(INSYM.EQ.1) THEN
  448.          ILUG=.TRUE.
  449.       ELSE
  450.          ILUG=.FALSE.
  451.       ENDIF
  452.       CALL ECCTAB(MTABLE,'MOT     ',4,XVA,'NSYM',ILOG,IOB,
  453.      $                   'LOGIQUE ',IOB,XVA,' ',ILUG,IOB)
  454.       if(idepe.ne.0) then
  455. *  on passe les ideme* a mrem sous forme de listchpo
  456.         n1=if
  457.         segini mlchpo,mlchp1
  458.         do i=1,if
  459.           mlchpo.ichpoi(i)=ideme0(1,i)
  460.           mlchp1.ichpoi(i)=ideme1(1,i)
  461.         enddo
  462.         CALL ECCTAB(MTABLE,'ENTIER  ',10,XVA,' ',ILOG,IOB,
  463.      $                   'LISTCHPO',IOB,XVA,' ',ILOG,mlchpo)
  464.         CALL ECCTAB(MTABLE,'ENTIER  ',11,XVA,' ',ILOG,IOB,
  465.      $                   'LISTCHPO',IOB,XVA,' ',ILOG,mlchp1)
  466. *  pour mrem on met la derniere raideur condensee. Elle contient les pointeurs pour remonter
  467.         CALL ECCTAB(MTABLE,'ENTIER  ',50,XVA,' ',ILOG,IOB,
  468.      $                   'RIGIDITE',IOB,XVA,' ',ILOG,ipoiri)
  469.       endif
  470.       SEGDES MRIGID
  471.       DO 26 I=IDEMEM(/1),1,-1
  472.          ISECO=IDEMEM(I)
  473.          CALL ACTOBJ ('CHPOINT ',ISECO,1)
  474.          CALL ECROBJ ('CHPOINT ',ISECO)
  475.   26  CONTINUE
  476.       SEGSUP IDEMEM
  477.       CALL ECROBJ ('TABLE   ',ISUPLO)
  478.       SEGINI MTEXTE
  479.       LTT=8
  480.       MTEXT(1:LTT) ='UNILATER'
  481.       NCART=8
  482.       SEGDES MTEXTE
  483.       CALL ECROBJ('TEXTE',MTEXTE)
  484.       mrigid=iprigo
  485.       segdes mrigid
  486.             norinc=norins
  487.       RETURN
  488.  
  489. C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ...
  490.    30 CONTINUE
  491. * il se peut que le dernier chp soit du frottement
  492. * on l'enleve car il ne sert a rien si on n'appele pas unilater
  493.         if (idemem(/1).gt.1.and.idepe.ne.0) then
  494.          mchpoi=ideme0(idemem(/1),if)
  495.          segact MCHPOI
  496.          if (mtypoi.eq.'LX  ') idemem(/1)=idemem(/1)-1
  497.         endif
  498. *  frottement
  499.       SEGDES IDEMEM
  500. *      write(6,*) ' ichski, igradj,insym ',ichski, igradj,insym
  501. *     write (6,*) ' imtvid ',imtvid
  502.       if (imtvid.eq.1) then
  503. *  matrice vide
  504. ***     write(6,*) ' attention matrice vide. Système surcontraint '
  505.         call erreur(-364)
  506. *
  507.         nsoupo=0
  508.         nat=0
  509.         segact idemem*mod
  510.  
  511.         do i=1,idemem(/1)
  512.          segini mchpoi
  513.          mchpoi.ifopoi = ifochs
  514.          idemem(i)=mchpoi
  515.         enddo
  516.         if (noen.eq.0) then
  517.           nat=2
  518.           nsoupo=0
  519.           segini mchpo4
  520.           call ecrobj('CHPOINT ',mchpo4)
  521.           call ecrent(0)
  522.         endif
  523.       else
  524. *       write(6,*) ' appel resou1 --  idemem(1)'
  525. *       segact idemem
  526. *       idesec= idemem(1)
  527. *       call ecchpo(idesec,0)
  528. *       write(6,*) ' appel resou1 -- ipoiri'
  529. *       call prrigi ( ipoiri,1)
  530. *       write(6,*) ' ichski  insym  ', ichski, insym
  531.        IF(ICHSKI.EQ.1) CALL RESOU1(IPOIRI,IDEMEM,NOID,NOEN,PREC,
  532.      >        ISTAB,ISOUCI,lagdua)
  533.        IF(IGRADJ.EQ.1) CALL GRACO0(IPOIRI,IDEMEM,NOID,NOEN)
  534.        IF(INSYM .EQ.1) CALL LDMT  (IPOIRI,IDEMEM,NOID,NOEN,PREC,ISOUCI,
  535.      >        LAGDUA )
  536.        IF(IERR.NE.0) GO TO 5001
  537.       endif
  538. C
  539. C-------------------------------------------------------
  540. C     LA SOLUTION EST CALCULEE --> ON LA MET EN FORME
  541.  
  542. **     if (noen.eq.0) then
  543. **       call lirobj('MAILLAGE',ipt8,1,iretou)
  544. **       if (ierr.ne.0) then
  545. **          norinc=norins
  546. **          return
  547. **       endif
  548. **       segact ipt8
  549. **       call lirent(nben,1,iretou)
  550. **       if (ierr.ne.0) return
  551. **     endif
  552.  
  553.       SEGACT IDEMEM*mod
  554.       N=IDEMEM(/1)
  555.       do i=1,n
  556.         mchpoi=idemem(i)
  557. *     les champs de points qui sortent sont de nature diffuse
  558.          SEGACT MCHPOI*MOD
  559.          ifopoi = ifochs
  560.          NAT = MAX(1,JATTRI(/1))
  561.          NSOUPO=IPCHP(/1)
  562.          SEGADJ MCHPOI
  563.          JATTRI(1)=1
  564.       enddo
  565.  
  566.       do 2010 ifois=1,30
  567.       segact mrigid
  568.       mrigid=jrsup
  569.       if (mrigid.eq.0) goto 2011
  570.       if(ierr.ne.0) then
  571.             norinc=norins
  572.             return
  573.          endif
  574.       call resour(idemem,ideme0,ideme1,mrigid,if,noen,ipt8,
  575.      >       isouci,iverif)
  576.       if (ierr.ne.0) then
  577.             norinc=norins
  578.             return
  579.          endif
  580.       if=if-1
  581.  2010 continue
  582.  2011 continue
  583. *
  584. *  on n'appelle plus verlx car je ne vois pas pourquoi on voudrait que les multiplicateurs de lagrange non éliminés soient nuls
  585. *
  586. ****   call verlx(ri2,iret,mchpo1,noen,ipt8)
  587. **     if (noen.eq.0) then
  588. **       call ecrent(nben)
  589. **     endif
  590. *
  591.       do 3 i=1,n
  592.        iret=idemem(n+1-i)
  593.  
  594. *     cas table de liaisons statiques
  595.       if (itbas.ne.0) then
  596.            il = n + 1 - i
  597.            ilo = idnote(il)
  598.          CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0,
  599.      &                   'TABLE',I1,X1,CHARRE,.true.,ITMOD)
  600.            if (ierr.ne.0) then
  601.             norinc=norins
  602.             return
  603.          endif
  604. c          CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
  605. c          CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
  606. c      &        'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
  607.  
  608.          CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE',
  609.      &        .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,IRET)
  610.  
  611.       else if (ipshpo.gt.0)  then
  612.          mlchpo = ipshpo
  613.          ichpoi(N+1-I) = iret
  614.       else
  615.          CALL ACTOBJ ('CHPOINT ',IRET,1)
  616.          CALL ECROBJ ('CHPOINT ',IRET)
  617.       endif
  618.  
  619.    3  CONTINUE
  620. c-----fin de boucle sur les solutions
  621.  
  622.  
  623. C-------------------------------------------------------
  624. c     MENAGE AVANT DE QUITTER
  625.  
  626.  5001 CONTINUE
  627.       if (itbas.ne.0) then
  628.          segdes mtab1
  629.          segsup idnote
  630.          CALL ECROBJ ('TABLE   ',itbas)
  631.       endif
  632.  
  633.       if (ipshpo.gt.0) then
  634.         mlchpo = ipshpo
  635.         CALL ACTOBJ ('LISTCHPO ',ipshpo,1)
  636.         CALL ECROBJ ('LISTCHPO ',ipshpo)
  637.       endif
  638.       SEGSUP IDEMEM
  639. C
  640.  5000 CONTINUE
  641.       norinc=norins
  642.       END
  643.  
  644.  
  645.  
  646.  
  647.  
  648.  
  649.  
  650.  
  651.  
  652.  

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