C CMCT6 SOURCE PV090527 26/04/30 21:15:20 12529 SUBROUTINE CMCT6(IRI1,IRI3) *_______________________________________________________________________ c c opérateur KPRE (cmct option ELEM) c c Réduction de la rigidité IRI1 sur la rigidité IRI3 C C Semblable a CMCT5 mais un peu plus general sans l'etre totalement C cf. Remarque plus loin c *_______________________________________________________________________ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) * -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMRIGID -INC SMELEME -INC CCHAMP -INC SMCOORD -INC SMCHPOI segment icpr(nbpts) segment inode(ino) segment jelnum(imaxel,ino) segment izone(imaxel,ino) segment INFMEL integer jmel,jcpr,jnode,kelnum,jzone endsegment * segment infzon integer iel1(nbsous) integer ielzo(nbsous) integer ilmasq(nele1,nbsous) integer jddlp(nbsous) integer jddld(nbsous) endsegment * segment ilmasq(nele1) segment iddlp(nligrp) segment iddld(nligrd) segment INFRIG INTEGER jzon(NRIGEL) * INTEGER jlmasq(NRIGEL) * INTEGER jddlp(NRIGEL) * INTEGER jddld(NRIGEL) endsegment logical ldescr,ldbg ldbg=.false. * * Le IMEL sur lequel on doit se baser est le maillage des points * atteints de IRI3, sachant que RI3 est deja partitionée comme le * maillage sous-tendant IRIGC2 * RI3=IRI3 * call rdscr2(ri3,ri4) IF (IERR.NE.0) RETURN NBNN=0 NBELEM=0 NBSOUS=RI4.IRIGEL(/2) NBREF=0 SEGINI MELEME do isous=1,nbsous lisous(isous)=RI4.IRIGEL(1,isous) enddo * Verif que le maillage ainsi constitue n'a pas d'elements en double ipt1=meleme call uniqma(ipt1,nbdif,0) IF (IERR.NE.0) RETURN if (ipt1.ne.meleme) segsup ipt1 if (nbdif.ne.0) then WRITE(IOIMP,*) $ 'CMCT6 Le maillage sous tendant RI3' $ ,' a des elements en double' call erreur(5) return endif segsup ri4 nbsou=lisous(/1) * * creation d'une numerotation locale * segini icpr ino=0 do 1 i =1, nbsou ipt4=lisous(i) segact ipt4 * write(ioimp,*) 'i=',i,' ipt4=',ipt4 * segprt,ipt4 do 2 j=1,ipt4.num(/2) do 22 k=1,ipt4.num(/1) ia= ipt4.num(k,j) if(icpr(ia).eq.0) then ino=ino+1 icpr(ia)=ino endif 22 continue 2 continue 1 continue * write(ioimp,*) 'ino=',ino * segprt,icpr * * on compte combien d'elements touche un noeud * segini inode do 3 i=1,nbsou ipt4=lisous(i) do 4 j=1,ipt4.num(/2) do 42 k=1,ipt4.num(/1) ia=ipt4.num(k,j) ib= icpr(ia) * write(ioimp,*) 'ia,ib=',ia,ib inode(ib)=inode(ib)+1 42 continue 4 continue 3 continue imaxel=0 do i=1,ino imaxel=max(imaxel,inode(i)) inode(i)=0 enddo segini jelnum,izone * * jelnum(i,k) dira le Ieme element qui touche le noeud k * izone (i,k) dira dans quel lisous le ieme element se trouve * attention le noeud K est le plus petit de l'élément * do 5 i=1,nbsou ipt4=lisous(i) do 6 j=1,ipt4.num(/2) ipeti=igrand do k=1,ipt4.num(/1) ipeti=min(ipt4.num(k,j),ipeti) enddo * write(ioimp,*) 'ipeti=',ipeti ib= icpr(ipeti) inode(ib)=inode(ib)+1 ic=inode(ib) jelnum(ic,ib)=j izone(ic,ib)=i 6 continue 5 continue * segprt,izone * segprt,jelnum * * Matrice réduite (sortie) * RI1=IRI1 SEGACT MELEME IF (IRI3.EQ.0) THEN write(ioimp,*) 'CMCT6 IRI3=0 non autorise' call erreur(5) return ENDIF RI3=IRI3 SEGACT RI3*MOD NRIG3=RI3.COERIG(/1) IF (NRIG3.NE.LISOUS(/1)) THEN write(ioimp,*) 'Erreur grave 1 dans cmct6' call erreur(5) return endif DO irig3=1,nrig3 DES3=RI3.IRIGEL(3,irig3) SEGACT,DES3 ENDDO RI3.MTYMAT='CMCT6' RI3.IFORIG=IFOUR RI3.MCRCNF=MCOORD DO irig3=1,nrig3 RI3.COERIG(irig3)=1.D0 *non ! ipt3=LISOUS(irig3) *non ! RI3.IRIGEL(1,irig3)=ipt3 ENDDO * * Traitement de la matrice d'entree * SEGACT RI1 CALL RDSCR2(RI1,RI2) IF (IERR.NE.0) RETURN RI1=RI2 NRIG1=RI1.COERIG(/1) NRIGEL=NRIG1 SEGINI INFRIG * DO IRIG1=1,nrig1 IPT1=RI1.IRIGEL(1,IRIG1) * Y a-t-il des elements a considerer SEGACT IPT1 * write(ioimp,*) 'IRIG1=',irig1, ' ipt1=',ipt1 * segprt,ipt1 DES1=RI1.IRIGEL(3,IRIG1) SEGACT DES1 nele1=ipt1.num(/2) nbsous=nbsou * segini ilmasq segini infzon nnod1=ipt1.num(/1) DO 11 iele1=1,nele1 ipeti=IGRAND DO inod1=1,nnod1 ipeti=min(ipeti,ipt1.num(inod1,iele1)) enddo * on regarde s'il existe un element de ipt1 ayant ce noeud en plus * petite position sinon on passe a l'élément suivant ib=icpr(ipeti) * write(ioimp,*) 'iele1=',iele1,' ipeti=',ipeti,' ib=',ib if(ib.eq.0) go to 11 if(inode(ib).eq.0) go to 11 do 13 mm=1,inode(ib) kzon=izone(mm,ib) ipt4 = lisous(kzon) if (ipt4.eq.ipt1) then ielzo(kzon)=ipt4.num(/2) ielzo(kzon)=ielzo(kzon)*(-1) if (iel1(kzon).eq.0) iel1(kzon)=iele1 goto 14 endif if (ipt4.num(/1).ne.ipt1.num(/1)) go to 13 iel=jelnum(mm,ib) do in=1,ipt1.num(/1) if (ipt1.num(in,iele1).ne.ipt4.num(in,iel))go to 13 enddo * on a trouve un element a conserver ielzo(kzon)=ielzo(kzon)+1 ilmasq(iele1,kzon)=iel if (iel1(kzon).eq.0) iel1(kzon)=iele1 goto 11 13 continue 11 continue 14 continue jzon(irig1)=infzon if (ldbg) then write(ioimp,*) 'irig1=',irig1,' nele1,nnod1=',nele1,nnod1 do izon=1,nbsous if (iel1(izon).ne.0) then il1=iel1(izon) il3=ilmasq(il1,izon) write(ioimp,*) ' izon=',izon,' ielzo=',ielzo(izon) $ ,' il1=',il1,' il3=',il3 else write(ioimp,*) ' izon=',izon,' ielzo=',ielzo(izon) endif enddo endif * segprt,RI3 do 21 izon=1,nbsous * tableau de correspondance des DESCR * write(ioimp,*) ' izon=',izon if (ielzo(izon).eq.0) goto 21 il1=iel1(izon) if (ielzo(izon).gt.0) then il3=ilmasq(il1,izon) else il3=il1 endif DES3=RI3.IRIGEL(3,IZON) ipt3=RI3.IRIGEL(1,izon) segact ipt3 DES1=RI1.IRIGEL(3,IRIG1) SEGACT DES1 * pour la primale NLIGP1=DES1.NOELEP(/1) NLIGRP=NLIGP1 SEGINI IDDLP DO ILIGP1=1,NLIGP1 * Remarque : * Bien foireux, on se base sur un seul element pour construire le * descr, il y a possibilite d'incoherence, il serait mieux d'etendre * la matrice RI1 sur les melemes de RI3 nno1=ipt1.num(des1.noelep(iligp1),il1) * write(ioimp,*) 'iligp1,nno1=',iligp1,nno1 ilig3=0 nligp3=des3.noelep(/1) do 113 iligp3=1,nligp3 nno2=ipt3.num(des3.noelep(iligp3),il3) * write(ioimp,*) 'iligp3,nno2=',iligp3,nno2 if (nno2.ne.nno1) goto 113 if (des3.lisinc(iligp3).ne.des1.lisinc(iligp1)) $ goto 113 ilig3=iligp3 goto 114 113 continue 114 continue IF (ilig3.EQ.0) THEN write(ioimp,*) 'Inconnue primale de nom : ' $ ,des1.lisinc(iligp1),' noeud : ' $ ,des1.noelep(iligp1) write(ioimp,*) ' presente dans la rigidite : ',ri1 write(ioimp,*) ' non presente dans la rigidite : ' $ ,ri3 write(ioimp,*) 'ipt1' segprt,ipt1 write(ioimp,*) 'des1' segprt,des1 write(ioimp,*) 'ipt3' segprt,ipt3 write(ioimp,*) 'des3' segprt,des3 call erreur(21) return else IDDLP(iligp1)=ILIG3 endif ENDDO JDDLP(izon)=IDDLP * segprt,iddlp * pour la duale NLIGD1=DES1.NOELED(/1) NLIGRD=NLIGD1 SEGINI IDDLD DO ILIGD1=1,NLIGD1 nno1=ipt1.num(des1.noeled(iligd1),il1) ilig3=0 nligd3=des3.noeled(/1) do 213 iligd3=1,nligd3 nno2=ipt3.num(des3.noeled(iligd3),il3) if (nno2.ne.nno1) goto 213 if (des3.lisdua(iligd3).ne.des1.lisdua(iligd1)) $ goto 213 ilig3=iligd3 goto 214 213 continue 214 continue IF (ilig3.EQ.0) THEN write(ioimp,*) 'Inconnue duale de nom : ' $ ,des1.lisdua(iligd1),' noeud : ' $ ,des1.noeled(iligd1) write(ioimp,*) ' presente dans la rigidite : ',ri1 write(ioimp,*) ' non presente dans la rigidite : ' $ ,ri3 write(ioimp,*) 'des1' segprt,des1 write(ioimp,*) 'des3' segprt,des3 call erreur(21) return else IDDLD(iligd1)=ILIG3 endif ENDDO JDDLD(izon)=IDDLD * segprt,iddld 21 continue ENDDO * do irig3=1,nrig3 * des3=ri3.irigel(3,irig3) * write(ioimp,*) 'irig3=',irig3 * segprt,des3 * enddo * * Transfert des xmatr selon le descripteur final des3 * DO IRIG1=1,nrig1 infzon=jzon(irig1) xc1=RI1.COERIG(IRIG1) * if (irig1.EQ.1) write(ioimp,*) 'xc1=',xc1 XMATR1=RI1.IRIGEL(4,IRIG1) SEGACT XMATR1 do 31 izon=1,nbsous if (ielzo(izon).eq.0) goto 31 * if (irig1.EQ.1) write(ioimp,*) 'izon=',izon iddlp=jddlp(izon) * if (irig1.EQ.1) segprt,iddlp iddld=jddld(izon) * if (irig1.EQ.1) segprt,iddld xmatr3=ri3.irigel(4,izon) if (xmatr3.eq.0) then ipt3=ri3.irigel(1,izon) des3=ri3.irigel(3,izon) NLIGRD=des3.NOELED(/1) NLIGRP=des3.NOELEP(/1) NELRIG=ipt3.num(/2) rigrel=0 SEGINI XMATR3 ri3.irigel(4,izon)=xmatr3 ri3.irigel(7,izon)=max(ri3.irigel(7,izon),ri1.irigel(7 $ ,irig1)) xmatr3.symre=max(xmatr3.symre,xmatr1.symre) endif NLIGD1=XMATR1.RE(/1) NLIGP1=XMATR1.RE(/2) NELE1 =XMATR1.RE(/3) if (ielzo(izon).gt.0) then DO iele1=1,nele1 iele3=ilmasq(iele1,izon) IF (iele3.ne.0) THEN DO ILIGP1=1,NLIGP1 ILIGP3=IDDLP(ILIGP1) DO ILIGD1=1,NLIGD1 ILIGD3=IDDLD(ILIGD1) XMATR3.RE(ILIGD3,ILIGP3,IELE3)= $ XMATR3.RE(ILIGD3,ILIGP3,IELE3) $ +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1 ENDDO ENDDO ENDIF ENDDO else DO iele1=1,nele1 DO ILIGP1=1,NLIGP1 ILIGP3=IDDLP(ILIGP1) DO ILIGD1=1,NLIGD1 ILIGD3=IDDLD(ILIGD1) XMATR3.RE(ILIGD3,ILIGP3,IELE1)= $ XMATR3.RE(ILIGD3,ILIGP3,IELE1) $ +XMATR1.RE(ILIGD1,ILIGP1,IELE1)*xc1 ENDDO ENDDO ENDDO endif 31 continue ENDDO * do irig3=1,nrig3 * xmatr3=ri3.irigel(4,irig3) * write(ioimp,*) 'irig3=',irig3 * if (xmatr3.gt.0) then * segprt,xmatr3 * else * write(ioimp,*) 'xmatr3=',xmatr3 * endif * enddo * Menage DO IRIG1=1,nrig1 infzon=jzon(irig1) do 41 izon=1,nbsous if (ielzo(izon).eq.0) goto 41 iddlp=jddlp(izon) if (iddlp.gt.0) segsup iddlp iddld=jddld(izon) if (iddld.gt.0) segsup iddld 41 continue segsup infzon ENDDO SEGSUP INFRIG segsup meleme segsup jelnum,izone segsup inode segsup icpr RETURN END