1 : C EXINCR SOURCE GOUNAND 11/06/14 21:15:22 7005 2 : SUBROUTINE EXINCR(RI1,MLMOT1,MLMOT2,RI2) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : EXINCR 7 : C DESCRIPTION : Extrait d'une rigidité la sous-matrice 8 : C d'inconnues primales et duales celles données 9 : C en argument CH*4 10 : * Dans le cas d'une relation, il serait prudent 11 : * de créer un nouveau multiplicateur mais alors si 12 : * on recombine les matrices extraites, on ne reobtient 13 : * pas l'originale. D'autre part, si on ne crée pas 14 : * de nouveau multiplicateur reso buggera sur la 15 : * matrice recombinée (présence d'un même multiplicateur 16 : * dans deux matrices différentes) On choisit de ne pas 17 : * créer de multiplicateur ici. 18 : C 19 : C 20 : C LANGAGE : ESOPE 21 : C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF) 22 : C mél : gounand@semt2.smts.cea.fr 23 : C*********************************************************************** 24 : C APPELES : 25 : C APPELES (E/S) : 26 : C APPELES (BLAS) : 27 : C APPELES (CALCUL) : 28 : C APPELE PAR : 29 : C*********************************************************************** 30 : C SYNTAXE GIBIANE : 31 : C RIG2 = EXTR RIG1 LMOT1 LMOT2 ; 32 : C 33 : C ENTREES : 34 : C ENTREES/SORTIES : 35 : C SORTIES : 36 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 37 : C*********************************************************************** 38 : C VERSION : v1, 01/06/2011, version initiale 39 : C HISTORIQUE : v1, 01/06/2011, création 40 : C HISTORIQUE : 41 : C HISTORIQUE : 42 : C*********************************************************************** 43 : -INC CCOPTIO 44 : -INC SMRIGID 45 : -INC SMELEME 46 : -INC SMLMOTS 47 : * 48 : * Segment de correspondance : noms d'inconnues 49 : * 50 : segment trav1 51 : integer cor1p(nligp2),cor1d(nligd2) 52 : endsegment 53 : segment trav2 54 : integer cor1n(nbnn2),cor2n(nbnn) 55 : endsegment 56 : * 57 : * Segment extensibles pour ranger les rigidités 58 : * créés 59 : * 60 : segment imatno(0) 61 : segment idesno(0) 62 : segment igeono(0) 63 : segment coerno(0) 64 : segment irg8no(0) 65 : segment irg7no(0) 66 : segment irg6no(0) 67 : segment irg5no(0) 68 : segment irg2no(0) 69 : * 70 : logical lsym,lrela 71 : * 72 : * Executable statements 73 : * 74 : SEGACT ri1 75 : SEGACT MLMOT1,MLMOT2 76 : * 77 : segini imatno,idesno,igeono,coerno,irg7no,irg5no 78 : segini irg2no,irg6no,irg8no 79 : do 1000 irig=1,ri1.irigel(/2) 80 : des1=ri1.irigel(3,irig) 81 : segact des1 82 : nligrp=des1.lisinc(/2) 83 : nligrd=des1.lisdua(/2) 84 : nligp2=nligrp 85 : nligd2=nligrd 86 : segini trav1 87 : * On regarde le descripteur et les inconnues concernées 88 : nligp2=0 89 : do 410 iligrp=1,nligrp 90 : do 420 incp=1,mlmot1.mots(/2) 91 : if (mlmot1.mots(incp).eq.des1.lisinc(iligrp)) then 92 : nligp2=nligp2+1 93 : cor1p(nligp2)=iligrp 94 : goto 410 95 : endif 96 : 420 continue 97 : 410 continue 98 : if (nligp2.eq.0) goto 999 99 : nligd2=0 100 : do 510 iligrd=1,nligrd 101 : do 520 incd=1,mlmot2.mots(/2) 102 : if (mlmot2.mots(incd).eq.des1.lisdua(iligrd)) then 103 : nligd2=nligd2+1 104 : cor1d(nligd2)=iligrd 105 : goto 510 106 : endif 107 : 520 continue 108 : 510 continue 109 : * WRITE(IOIMP,*) 'irig=',irig 110 : * WRITE(IOIMP,*) 'cor1p=' 111 : * WRITE (IOIMP,2020) (cor1p(I),I=1,cor1p(/1)) 112 : * WRITE(IOIMP,*) 'cor1d=' 113 : * WRITE (IOIMP,2020) (cor1d(I),I=1,cor1d(/1)) 114 : * 2020 FORMAT (20(2X,I4) ) 115 : if (nligd2.eq.0) goto 999 116 : * On fait un cas particulier pour les matrices de relations 117 : * Il faut vérifier qu'on a bien gardé la symétrie. 118 : * Il faut également que la relation ait gardé le multiplicateur 119 : * et porte au moins sur un ddl autre ! 120 : ipt1=ri1.irigel(1,irig) 121 : segact ipt1 122 : lrela=.false. 123 : if (ipt1.itypel.eq.22) then 124 : if (nligp2.eq.nligd2) then 125 : do 700 ilig=1,nligp2 126 : if (cor1p(ilig).ne.cor1d(ilig)) goto 710 127 : 700 continue 128 : lrela=.true. 129 : 710 continue 130 : endif 131 : endif 132 : * l'ordre de ces deux lignes est important 133 : if (lrela.and.cor1p(1).ne.1) goto 998 134 : if (lrela.and.nligp2.eq.1) goto 998 135 : coerno(**)=ri1.coerig(irig) 136 : irg2no(**) =ri1.irigel(2,irig) 137 : irg5no(**) =ri1.irigel(5,irig) 138 : irg6no(**) =ri1.irigel(6,irig) 139 : irg8no(**) =ri1.irigel(8,irig) 140 : * Le cas particulier facile, où on garde toute la matrice 141 : * élémentaire : on recopie la matrice originale 142 : if (nligp2.eq.nligrp.and.nligd2.eq.nligrd) then 143 : igeono(**)=ri1.irigel(1,irig) 144 : idesno(**) =des1 145 : imatno(**) =ri1.irigel(4,irig) 146 : irg7no(**) =ri1.irigel(7,irig) 147 : goto 999 148 : endif 149 : * Sinon, il faut réduire la matrice et son descripteur, voire même 150 : * la géométrie 151 : * Réduction du descripteur 152 : nligrp=nligp2 153 : nligrd=nligd2 154 : segini des2 155 : do iligp2=1,nligp2 156 : iligp1=cor1p(iligp2) 157 : des2.lisinc(iligp2)=des1.lisinc(iligp1) 158 : des2.noelep(iligp2)=des1.noelep(iligp1) 159 : enddo 160 : do iligd2=1,nligd2 161 : iligd1=cor1d(iligd2) 162 : des2.lisdua(iligd2)=des1.lisdua(iligd1) 163 : des2.noeled(iligd2)=des1.noeled(iligd1) 164 : enddo 165 : idesno(**)=des2 166 : * Réduction de la matrice 167 : xmatr1=ri1.irigel(4,irig) 168 : segact xmatr1 169 : nelrig=xmatr1.re(/3) 170 : segini xmatr2 171 : do ielrig=1,nelrig 172 : do iligp2=1,nligp2 173 : iligp1=cor1p(iligp2) 174 : do iligd2=1,nligd2 175 : iligd1=cor1d(iligd2) 176 : xmatr2.re(iligd2,iligp2,ielrig)= 177 : $ xmatr1.re(iligd1,iligp1,ielrig) 178 : enddo 179 : enddo 180 : enddo 181 : segdes xmatr2 182 : segdes xmatr1 183 : imatno(**) =xmatr2 184 : * Réduction éventuelle de la géométrie si des noeuds 185 : * ne sont pas référencés 186 : nbnn=ipt1.num(/1) 187 : nbnn2=ipt1.num(/1) 188 : segini trav2 189 : nbnn2=0 190 : do in=1,des2.noelep(/1) 191 : ic1=cor2n(des2.noelep(in)) 192 : if (ic1.eq.0) then 193 : nbnn2=nbnn2+1 194 : cor1n(nbnn2)=des2.noelep(in) 195 : cor2n(des2.noelep(in))=nbnn2 196 : endif 197 : enddo 198 : do in=1,des2.noeled(/1) 199 : ic1=cor2n(des2.noeled(in)) 200 : if (ic1.eq.0) then 201 : nbnn2=nbnn2+1 202 : cor1n(nbnn2)=des2.noeled(in) 203 : cor2n(des2.noeled(in))=nbnn2 204 : endif 205 : enddo 206 : * Si tous les noeuds ont été vus, on peut peut-être 207 : * garder la géométrie 208 : if (nbnn2.eq.nbnn) then 209 : if (ipt1.itypel.eq.22.and.(.not.lrela)) then 210 : segini,ipt2=ipt1 211 : ipt2.itypel=28 212 : segdes ipt2 213 : else 214 : ipt2=ipt1 215 : endif 216 : else 217 : * Sinon, on ne garde que les noeuds vus dans le nouveau 218 : * maillage et on modifie le descripteur 219 : nbnn=nbnn2 220 : nbelem=ipt1.num(/1) 221 : nbsous=0 222 : nbref=0 223 : segini ipt2 224 : ipt2.itypel=28 225 : if (lrela) ipt2.itypel=22 226 : do iel=1,nbelem 227 : do ibnn2=1,nbnn2 228 : ibnn=cor1n(ibnn2) 229 : ipt2.num(ibnn2,iel)=ipt1.num(ibnn,iel) 230 : enddo 231 : enddo 232 : do iligrp=1,des2.noelep(/1) 233 : des2.noelep(iligrp)=cor2n(des2.noelep(iligrp)) 234 : enddo 235 : do iligrd=1,des2.noeled(/1) 236 : des2.noeled(iligrd)=cor2n(des2.noeled(iligrd)) 237 : enddo 238 : segdes ipt2 239 : endif 240 : segsup trav2 241 : igeono(**)=ipt2 242 : segdes des2 243 : * Garde-t-on la symétrie ou l'antisymétrie ? Par défaut non, mais si 244 : * les ddl primaux et duaux que l'on garde sont les mêmes, oui ! 245 : lsym=.false. 246 : if (ri1.irigel(7,irig).ne.2) then 247 : if (nligp2.eq.nligd2) then 248 : do 600 ilig=1,nligp2 249 : if (cor1p(ilig).ne.cor1d(ilig)) goto 610 250 : 600 continue 251 : lsym=.true. 252 : 610 continue 253 : endif 254 : endif 255 : if (lsym) then 256 : irg7no(**) =ri1.irigel(7,irig) 257 : else 258 : irg7no(**) =2 259 : endif 260 : 998 continue 261 : segdes ipt1 262 : 999 continue 263 : segdes des1 264 : segsup trav1 265 : 1000 CONTINUE 266 : SEGDES MLMOT1,MLMOT2 267 : * 268 : * il ne reste plus qu'a ranger dans ri2 les nouvelles raideurs engendrées 269 : * 270 : NRIGEL=imatno(/1) 271 : SEGINI RI2 272 : RI2.MTYMAT=RI1.MTYMAT 273 : SEGDES RI1 274 : DO irig=1,NRIGEL 275 : RI2.COERIG(irig)=coerno(irig) 276 : RI2.IRIGEL(1,irig)=igeono(irig) 277 : RI2.IRIGEL(2,irig)=irg2no(irig) 278 : RI2.IRIGEL(3,irig)=idesno(irig) 279 : RI2.IRIGEL(4,irig)=imatno(irig) 280 : RI2.IRIGEL(5,irig)=irg5no(irig) 281 : RI2.IRIGEL(6,irig)=irg6no(irig) 282 : RI2.IRIGEL(7,irig)=irg7no(irig) 283 : RI2.IRIGEL(8,irig)=irg8no(irig) 284 : ENDDO 285 : SEGDES RI2 286 : segsup imatno,idesno,igeono,coerno,irg7no,irg5no 287 : segsup irg2no,irg6no,irg8no 288 : * 289 : * Normal termination 290 : * 291 : * IRET=0 292 : RETURN 293 : * 294 : * End of subroutine EXINCR 295 : * 296 : END 297 : 298 : 299 : 300 : 301 : 302 :
© Cast3M 2003 - All rights reserved.
Disclaimer