Sources Esope | Cast3M

kres6
C KRES6     SOURCE    PV090527  23/09/13    21:15:02     11739          
      SUBROUTINE KRES6(MRIGID,KSMBR,LDMULT,NELIM,
     $     MRIGIC,KSMBRC,KSMBR1)
*      SUBROUTINE KRES6(MRIGID,KSMBR,IDEPE,
*     $     MRIGIC,KSMBRC,KSMBR0,KSMBR1)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : KRES6
C DESCRIPTION : Effectue la condensation des relations
C               Repris de resou.eso
C
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C VERSION    : v1, 15/06/2011, version initiale
C HISTORIQUE : v1, 15/06/2011, création
C HISTORIQUE : 2019/04/10 remplacement de NOEL par NELIM
C     Idéalement, il faudrait reprendre ce que Pierre a fait dans
C     RESOU dans les fiches 10[0,1]?? et avec MREM.
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC SMRIGID
-INC SMCHPOI
-INC SMELEME
*
*  Logique indiquant si on dualise les multiplicateurs de Lagrange
      LOGICAL LDMULT
*
      NOUNIL=0
      NOEN=1
      IPOIRI=MRIGID
*  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
      nrige= irigel(/1)
      idepe=0
      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+1
         if (irigel(6,irig).ne.0) iunil=1
         segdes meleme
 1000 continue
*      idepe=0
      lagdua=0
      if(idepe.ne.0) then
C  on va separer les raideurs
*         write (6,*) ' nounil jrcond iunil idepe vaut ',nounil,jrcond,
*     $        iunil, idepe
         if (jrcond.eq.0) then
*transmis en argument            nelim=1
            if (ldmult) then
               call separm(mrigid,ri1,ri2,nounil,lagdua,nelim,0)
            else
               call separm(mrigid,ri1,ri2,nounil,-1,nelim,0)
            endif
*            if (lagdua.ne.0) then
*               write(6,*) ' resou apres separm -- lagdua'
*               call ecmail(lagdua)
*            else
*               write(6,*) ' resou apres separm -- lagdua=0'
*            endif
*            write(6,*) ' resou apres separm -- ri2'
*            call prrigi(ri2,1)
*            write(6,*) ' resou apres separm -- ri1'
*            call prrigi ( ri1,0)
*            write(6,*) ' fin impression ri1'
            segact mrigid*mod
            jrelim=ri1
            jrgard=ri2
            imlag=lagdua
            call fusrig(ri1,ri2,ipoir0)
            jrtot=ipoir0
         else
            ri1=jrelim
            ri2=jrgard
            ipoir0=jrtot
            lagdua=imlag
            ipt1=lagdua
            if (ipt1.ne.0) segact ipt1
         endif
C
*    mrigid matrice complete
*    ri1    dependance
*    ri2    les autres matrices
*    ri6    matrice de transfert
*    ri3    matrice reduite
*    ri5    matrice de transfert transposee
C
C   on va proceder  a la condensation rigidite puis forces
         if(jrcond.eq.0) then
            CALL DEPEN3(RI1,RI6)
*            write (6,*) ' resou apres depen3 --- ri6 '
*            call prrigi(ri6,1)
            call scnd2 (ri2,ri6,ri3)
*            write (6,*) '  '
*            write (6,*) '  '
*            write (6,*) ' resou  apres scnd2--- ri3 '
*            write (6,*) ' '
*            call prrigi(ri3,1)
            segact ri3
            if(ierr.ne.0) then
               segsup ri1,ri2,ri6
               return
            endif
 
            segact mrigid*mod
            jrcond=ri3
            JRDEPP=RI6
C     dualisation  de la (les) matrice(s) de dependance
            call dual00(ri6,ri5)
*            write(6,*) ' apres dual0  -- ri5'
*            call prrigi( ri5,1)
            jrdepd=ri5
            ipoiri = ri3
         else
            ipoiri= jrcond
            RI6   = JRDEPP
            ri5   = jrdepd
         endif
*  test si ri3 est vide
         ri3=jrcond
         segact ri3
*        write (6,*) ' dans resou ri3.irigel(/2) ',ri3.irigel(/2)
         if (ri3.irigel(/2).eq.0) imtvid=1
C
         segdes ri1,ri2,mrigid
         mrigic=ipoiri
C   maintenant  les seconds membres
C        write(6,*) ' ipoiri jrdepp jrdepd',ipoiri,ri6,ri5
C         call prrigi(ri3)
C        call prrigi(ri5)
* en cas de contacts on ne dualise pas . Ce sera fait dans unilater
         ifrot=0
         MRIGID=IPOIRI
         SEGACT MRIGID*MOD
         DO  I=1,IRIGEL(/2)
            IF(IRIGEL(6,I).ne.0) ifrot=1
         enddo
         if (nounil.eq.1) ifrot=0
*     if (ifrot.eq.0) write (6,*) ' resou non unilateral '
*     if (ifrot.eq.1) write (6,*) ' resou unilateral pas de dualisation'
         if (lagdua.ne.0) then
            ipt8=lagdua
            segact ipt8
*     call ecmail(lagdua,0)
         endif
*
         ichp2=ksmbr
*         ksmbr0=ichp2
*  transferer les valeurs imposees des relations sur les inconnues (à éliminer
C )
*            write (6,*) ' valeurs imposees ichp2'
*            call ecchpo(ichp2,0)
         call transr(ichp2,ri1,ichp3)
*            call prrigi(ri1)
*            write (6,*) ' apres transfert ichp3'
*            call ecchpo(ichp3,0)
         ksmbr1=ichp3
         call mucpri(ichp3,ri2,ichp4)
*            write (6,*) ' apres mucpri ichp4 '
*            call ecchpo(ichp4,0)
*  ri2 est deja dualise. Il faut donc dedualiser ichp4
         call dbbcd(ichp4,lagdua)
*            write (6,*) ' apres dbbcd ichp4 '
*            call ecchpo(ichp4,0)
         call adchpo(ichp2,ichp4,ichp1,1.D0,-1.0D0)
         call dtchpo(ichp4)
*            write (6,*) ' apres  adchpo ichp1'
*            call ecchpo(ichp1,0)
         call mucpri(ichp1,ri5,ichp2)
*            write (6,*) ' impression de ri5 '
*            call prrigi(ri5,1)
*            write (6,*) ' apres mucpri ichp2'
*            call ecchpo(ichp2,0)
C         mchpo1=ichp1
C         segact mchpo1
C         write(6,*) 'reso mchpo1 jattri(1)',mchpo1.jattri(1)
C         segdes mchpo1
C
         mchpo2= ichp1
         segact mchpo2*mod
         mchpo2.jattri(1)=2
         mchpo2= ichp2
         segact mchpo2*mod
         mchpo2.jattri(1)=2
C         write(6,*) 'reso mchpo2 jattri(1)',mchpo2.jattri(1)
         segdes mchpo2
C
         call fuchpo (ichp1,ichp2,ichp3)
*     call dtchpo(ichp1)
         call dtchpo(ichp2)
* Ajout gounand : à ce stade, la force réduite n'est pas nulle sur les
* ddls supprimés (multiplicateurs de Lagrange et ddl dépendants), on les
* enlève.
*  vecteur a resoudre
*            write (6,*) ' le vecteur avant reduction '
*            call ecchpo(ichp3,0)
         ichp2=ichp3
         CALL redfor(ichp2,ri1,ichp3)
         if (ierr.ne.0) return
         call dtchpo(ichp2)
         ksmbrc=ichp3
*  vecteur a resoudre
*            write (6,*) ' le vecteur '
*            call ecchpo(ichp3,0)
*  dualisation des mult de lagrange
         if (lagdua.ne.0.and.ifrot.eq.0) call dbbch(ichp3,lagdua)
*  matrice
*            write (6,*) ' la matrice '
*            call prrigi(ipoiri)
         noid = 1
      else
         mrigic=mrigid
         ksmbrc=ksmbr
         ksmbr1=0
      endif
 
 
 
      RETURN
*
* End of subroutine KRES6
*
      END