Esope source files

   1 : C RELABA    SOURCE    CHAT      09/10/09    21:22:39     6519           
   2 :       SUBROUTINE RELABA
   3 : *
   4 : * relation barycentrique pour les ddl d'un noeud vis-à-vis
   5 : * d'un maillage
   6 : *  rig1 = rela bary lmot1 (ou 'DEPL' ou 'ROTA') p1 geo1 
   7 : *
   8 : *
   9 :       IMPLICIT INTEGER(I-N)
  10 :       IMPLICIT REAL*8 (A-H,O-Z)
  11 : -INC CCOPTIO
  12 : -INC CCGEOME
  13 : -INC SMELEME
  14 : -INC SMCOORD
  15 : -INC SMRIGID
  16 : -INC CCHAMP
  17 : -INC SMLREEL
  18 : -INC SMLMOTS
  19 : c
  20 :       CHARACTER*4 ITCORI(2),LEMOT,LESDDL(6),LESDUA(6)
  21 :       CHARACTER*4 MODEPL(5),MODEDU(5),MOROTA(5),MORODU(5)
  22 :       DATA ITCORI(1)/'DEPL'/
  23 :       DATA ITCORI(2)/'ROTA'/
  24 :       DATA MODEPL / 'UX  ','UY  ','UZ  ','UR  ','UT  ' /
  25 :       DATA MODEDU / 'FX  ','FY  ','FZ  ','FR  ','FT  ' /
  26 :       DATA MOROTA / 'RX  ','RY  ','RZ  ','RT  ','RS  ' /
  27 :       DATA MORODU / 'MX  ','MY  ','MZ  ','MT  ','MS  ' /
  28 : 
  28 :       CALL LIROBJ('LISTMOTS',IP0,0,IRETOU)
  29 : 
  30 :       if (IP0.eq.0) then
  30 :          CALL LIRMOT(ITCORI,2,IMOT,0)
  31 :          if (imot.eq.0) then
  31 :           call erreur(5)
  32 :           return
  33 :          else if (imot.eq.1) then
  34 :           if (ifour.eq.2) then
  35 :            nosddl = 3 
  36 :            lesddl(1) = modepl(1)
  37 :            lesddl(2) = modepl(2)
  38 :            lesddl(3) = modepl(3)
  39 :           else
  40 :           endif
  40 :            CALL LIRMOT(ITCORI,2,IMOT,0)
  41 :            if (imot.eq.2) then
  42 :             if (ifour.eq.2) then
  43 :              nosddl = 6
  44 :              lesddl(4) = morota(1)
  45 :              lesddl(5) = morota(2)
  46 :              lesddl(6) = morota(3)
  47 :             else
  48 :             endif
  49 :            endif
  50 :          else
  50 :            CALL LIRMOT(ITCORI,2,IMOT,0)
  51 :            if (imot.eq.1) then
  52 :             if (ifour.eq.2) then
  53 :              nosddl = 6
  54 :              iddl = 3
  55 :             else
  56 :             endif
  57 :            else
  58 :             if (ifour.eq.2) then
  59 :              nosddl = 3
  60 :              iddl = 0
  61 :             else
  62 :             endif
  63 :             
  64 :            endif
  65 :            if (ifour.eq.2) then 
  66 :             lesddl(iddl + 1) = morota(1)
  67 :             lesddl(iddl + 2) = morota(2)
  68 :             lesddl(iddl + 3) = morota(3)
  69 :           else
  70 :           endif
  71 :         endif
  72 :       else
  73 :        mlmots = ip0
  74 :        segact mlmots
  75 :        jgn = mots(/2)
  76 :        do 11 imo = 1,jgn
  77 :         do jde =1,5
  78 :          if (mots(imo).eq.modepl(jde)) then
  79 :            lesddl(imo) = modepl(jde)
  80 :            lesdua(imo) = modedu(jde)
  81 :            goto 11
  82 :          endif
  83 :         enddo
  84 :         do jde = 1,5
  85 :          if (mots(imo).eq.morota(jde)) then
  86 :            lesddl(imo) = morota(jde)
  87 :            lesdua(imo) = morodu(jde)
  88 :            goto 11
  89 :          endif
  90 :         enddo
  91 : c le mot  n est pas un ddl 
  91 :         call erreur(5)
  92 :         return
  93 :  11    continue
  94 :        nosddl = jgn
  95 :       endif
  96 : 
  97 : 
  97 :       CALL LIROBJ('POINT',IP1,1,IRETOU)
  97 :       CALL LIROBJ('MAILLAGE',IP2,1,IRETOU)
  98 :       meleme = ip2
  99 :       segact meleme
  99 :       if (itypel.ne.1) call change(ip2,1)
 100 :       meleme = ip2
 101 :       segact meleme
 102 :       nbele1 = num(/2)
 103 :       ipt1 = ip2
 103 :       CALL LIROBJ('LISTREEL',IP3,0,IRETOU)
 104 :       if (iretou.eq.1) then
 105 :        mlreel = ip3      
 106 :        segact mlreel
 107 :        jg = prog(/1)
 108 :        if (jg.ne.nbele1) then
 109 : c autant de coef que de points
 109 :          call erreur(5)
 110 :          return
 111 :        endif
 112 :        som1 = prog(1)
 113 :        do j = 2,jg
 114 :          som1 = som1 + prog(j)
 115 :        enddo
 116 :        if (som1.eq.0.D0) then
 117 : c somme des coefs nulle
 117 :          call erreur(5)
 118 :          return
 119 :        endif 
 120 :       else
 121 :        jg = 1
 122 :       endif
 123 :       
 124 : 
 125 : * cree multiplicateur(s)
 126 :       NBPTI=XCOOR(/1)/(IDIM+1)
 127 :       NBPTS = NBPTI + nosddl
 128 :       segadj mcoord
 129 :       do j = 1,nosddl
 130 :        xcoor((nbpti - 1 + j)*(idim + 1) + 1) = 0.d0
 131 :        xcoor((nbpti - 1 + j)*(idim + 1) + 2) = 0.d0
 132 :        if (idim.eq.3) xcoor((nbpti - 1 + j)*(idim + 1) + 3) = 0.d0
 133 :        xcoor((nbpti + j)*(idim + 1)) = 0.d0
 134 :       enddo
 135 : 
 136 : * cree rigidite
 137 :       NRIGE = 6
 138 :       NRIGEL = nosddl
 139 :       segini mrigid
 140 :       ICHOLE=0
 141 :       IMGEO1=0
 142 :       IMGEO2=0
 143 :       ISUPEQ=0
 144 :       IFORIG=IFOMOD
 145 :       MTYMAT='RIGIDITE'
 146 : 
 147 :        NLIGRP=nbele1 + 2
 148 :        NLIGRD=nbele1 + 2
 149 :            nelrig=1
 150 :        segini xmatri
 151 : *       nelrig = 1
 152 : *       segini imatri
 153 : *       imattt(1) = xmatri
 154 :        re(1,1,1) = 0.d0
 155 :        if (jg.eq.1) then
 156 :          re(1,2,1) = 1.D0/nbele1
 157 :          re(2,1,1) = re(1,2,1)
 158 :          do l = 1, nbele1
 159 :           re(1,2 + l,1) = -1.D0
 160 :           re(2+l,1,1) = -1.D0
 161 :          enddo
 162 :        else
 163 :          re(1,2,1) = 1.D0/som1
 164 :          re(2,1,1) = re(1,2,1)
 165 :          do l = 1, nbele1
 166 :           re(1,2 + l,1) = -1.D0*prog(l)
 167 :           re(2+l,1,1) = re(1,2+l,1)
 168 :          enddo
 169 :        endif
 170 :        segdes xmatri
 171 : 
 172 :       do ir = 1,nosddl
 173 :         ipmu = NBPTI + ir
 174 :         nbnn = nbele1 + 2
 175 :         nbelem = 1
 176 :         nbsous = 0
 177 :         nbref = 0 
 178 :         segini meleme
 179 :         itypel = 22
 180 :         num(1,1) = ipmu
 181 :         num(2,1) = ip1
 182 :         segact ipt1
 183 :         do in = 1,nbele1
 184 :           num(2+in,1) = ipt1.num(1,in)
 185 :         enddo
 186 :         icolor(1) = idcoul
 187 :         segdes meleme
 188 : 
 189 :        SEGINI DESCR
 190 :        NOELEP(1)= 1
 191 :        NOELED(1)= 1
 192 :        LISINC(1)='LX'
 193 :        LISDUA(1)='FLX'
 194 :        do l = 1,nbele1+1
 195 :          LISINC(1 + l)=lesddl(ir)
 196 :          LISDUA(1 + l)=lesdua(ir)
 197 :          NOELEP(1+l)= 1 + l
 198 :          NOELED(1+l)= 1 + l
 199 :        enddo
 200 :        SEGDES DESCR
 201 : 
 202 :         coerig(ir) = 1.d0
 203 :         IRIGEL(2,ir) = 0
 204 :         IRIGEL(5,ir) = NIFOUR
 205 :         IRIGEL(6,ir) = 0
 206 :         IRIGEL(1,ir) = meleme
 207 :         IRIGEL(3,ir) = descr
 208 :         IRIGEL(4,ir) = xmatri
 209 : 
 210 : 
 211 :       enddo
 212 :       segdes mrigid
 213 : 
 213 :       CALL ECROBJ('RIGIDITE',mrigid)
 214 :       RETURN
 215 :       END
 216 :  
 217 :