cmct4
C CMCT4 SOURCE GOUNAND 25/06/11 21:15:04 12278 SUBROUTINE CMCT4(IRIGC,IRIGD,IRIGB,IRIGE,IMEL,IRIG2) *_______________________________________________________________________ c c opérateur KPRE (cmct option ELEM) c c entrée c IRIGB : rigidité B c IRIGC : rigidité C c IRIGD : rigidité D c IMEL : maillage c c sortie c IRIG2 : rigidité contenant la matrice C* D*-1 B*t construite par c elements sur les elements de MELEME c avec C*, D*, B* les matrices reduites sur les elements de MELEME c *_______________________________________________________________________ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) * -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMRIGID POINTEUR IRIGC.MRIGID,IRIGD.MRIGID,IRIGB.MRIGID POINTEUR IRIGC2.MRIGID,IRIGD2.MRIGID,IRIGB2.MRIGID POINTEUR IRIGC3.MRIGID,IRIGB3.MRIGID POINTEUR IRIGE.MRIGID,IRIGE2.MRIGID POINTEUR IRIG2.MRIGID POINTEUR XMATRC.XMATRI,XMATRD.XMATRI,XMATRB.XMATRI POINTEUR XMATRE.XMATRI POINTEUR XMAT.XMATRI POINTEUR XMATV.XMATRI SEGMENT vvec(ntot) SEGMENT amat(ntot,ntot) SEGMENT amats(ntot,ntot) SEGMENT amatv(ntot,ntot) -INC SMELEME -INC CCHAMP -INC SMCOORD -INC SMCHPOI * CHARACTER*(LOCHPO) NOMINC,NOMDUA * segment icpr(nbpts) segment inode(ino) segment jelnum(imaxel,ino) segment izone(imaxel,ino) segment INFMEL integer jmel,jcpr,jnode,kelnum,jzone endsegment logical lvdec,ldbg * ldbg=.false. if (ldbg )write(ioimp,*) 'irigc,irigd,irigb,irige=',irigc,irigd $ ,irigb,irige lvdec=.false. nprint=40 iprint=0 npivhg=0 npivbd=0 nelpnu=0 xzpiv=xzprec**0.75d0 xzfac=xzprec**0.5d0 IF (IMEL.NE.0) THEN if (ierr.ne.0) return irigc=irigc2 ENDIF * segact irigc IFOC=IRIGC.IFORIG IFO2=IFOC nbsou=irigc.irigel(/2) C Cas d'un MELEME vide => RIGIDITE vide IF (nbsou.eq.0) THEN NRIGEL = 0 SEGINI,IRIG2 IRIG2.MTYMAT = 'KPRE' IRIG2.IFORIG = IFO2 RETURN ENDIF C C On commence par régulariser irigc au sens ou tous les points de C son maillage sont atteints C if (ldbg) write(ioimp,*) 'IRIGC->MELEME' call rdscr2(irigc,irigc2) IF (IERR.NE.0) RETURN call melrig(irigc2,ipt5) IF (IERR.NE.0) RETURN IF (IERR.NE.0) RETURN C segact meleme nbsou=lisous(/1) SEGINI INFMEL INFMEL.JMEL=MELEME * * Réduction de C sur le maillage MELEME * if (ldbg) write(ioimp,*) ' IRIGC ' * call ecrobj('RIGIDITE',irigc) * call ecrcha('RESU') * call prlist IRIGC3=0 CALL CMCT5(IRIGC2,INFMEL,IRIGC3) IF (IERR.NE.0) RETURN segsup irigc2 irigc2=irigc3 * write(ioimp,*) ' IRIGC2 ' * segprt,irigc2 * call ecrobj('RIGIDITE',irigc2) * call ecrcha('RESU') * call prlist * return * * Réduction eventuelle de B sur le maillage MELEME et sur le DESCR de C2 * IF (IRIGB.NE.IRIGC) THEN SEGACT IRIGB IFOB=IRIGB.IFORIG IF (IFOB.NE.IFOC) THEN moterr(1:8)='RIGIDITE' interr(1)=IFOB interr(2)=IFOC interr(3)=IFOUR IFO2 = IFOUR ENDIF call rdscr2(irigb,irigb2) IF (IERR.NE.0) RETURN * Amélioration possible : transposition auto de B si necessaire if (ldbg) write(ioimp,*) ' IRIGB .NE. IRIGC' segini,irigb3=irigc2 * On enleve les xmatri nrig=irigb3.coerig(/1) do irig=1,nrig irigb3.irigel(4,irig)=0 enddo CALL CMCT5(IRIGB2,INFMEL,IRIGB3) IF (IERR.NE.0) RETURN segsup irigb2 irigb2=irigb3 * write(ioimp,*) ' IRIGB2 ' * segprt,irigb2 * call ecrobj('RIGIDITE',irigb2) * call ecrcha('RESU') * call prlist ELSE if (ldbg) write(ioimp,*) ' IRIGB = IRIGC' ENDIF * * Reduction de D sur le maillage MELEME et sur le DESCR primal de C2 * SEGACT IRIGD IFOD=IRIGD.IFORIG IF (IFOD.NE.IFOC) THEN moterr(1:8)='RIGIDITE' interr(1)=IFOD interr(2)=IFOC interr(3)=IFOUR IFO2 = IFOUR ENDIF NRIGEL=nbsou SEGINI,IRIGD2 SEGACT,IRIGC2 DO irig=1,NRIGEL ipt1=irigc2.irigel(1,irig) des1=irigc2.irigel(3,irig) call descar(des1,1,des3) IF (IERR.NE.0) RETURN irigd2.irigel(1,irig)=ipt1 irigd2.irigel(3,irig)=des3 ENDDO if (ldbg) write(ioimp,*) ' IRIGD ' CALL CMCT6(IRIGD,IRIGD2) IF (IERR.NE.0) RETURN * call ecrobj('RIGIDITE',irigd2) * call ecrcha('RESU') * call prlist * return * * Reduction (eventuelle) de E sur le maillage MELEME et sur le DESCR dual de C2 * IF (IRIGE.NE.0) THEN SEGACT IRIGD IFOE=IRIGE.IFORIG IF (IFOE.NE.IFOC) THEN moterr(1:8)='RIGIDITE' interr(1)=IFOE interr(2)=IFOC interr(3)=IFOUR IFO2 = IFOUR ENDIF NRIGEL=nbsou SEGINI,IRIGE2 SEGACT,IRIGC2 DO irig=1,NRIGEL ipt1=irigc2.irigel(1,irig) des1=irigc2.irigel(3,irig) call descar(des1,2,des3) IF (IERR.NE.0) RETURN irige2.irigel(1,irig)=ipt1 irige2.irigel(3,irig)=des3 ENDDO if (ldbg) write(ioimp,*) ' IRIGE ' CALL CMCT6(IRIGE,IRIGE2) IF (IERR.NE.0) RETURN ELSE IRIGE2=0 ENDIF * * On est pret maintenant à calculer E - C D-1 Bt * On envisage trois algorithmes : * 1) Force Brute : inversion de D puis produit * 2) Factorisation de la matrice (D Bt) = (Lu 0 ) ( Uu Upu) * (C 0 ) (Lup Lp) ( 0 Up ) * Si je comprends bien, Lp * Up = C D-1 Bt Oui ! * 3) Brute brute : d'après la formulation d'inversion, on inverse * tout, puis on inverse le bloc pp de l'inverse. Ca nous redonne le * complément de Schur * * * Algorithme factorisation LDMt: on pourrait pivoter symétriquement ? * Repris de Golub et Van Loan ! * NRIGEL=nbsou SEGINI IRIG2 IRIG2.MTYMAT = 'KPRE' IRIG2.IFORIG = IFO2 SEGACT IRIGC2 DO irig=1,NRIGEL IRIG2.COERIG(irig)=1.D0 ipt1=irigc2.irigel(1,irig) des1=irigc2.irigel(3,irig) nligd=des1.noeled(/1) nligp=des1.noelep(/1) call descar(des1,2,des3) IF (IERR.NE.0) RETURN nlrig=ipt1.num(/2) nelpnu=nelpnu+nlrig * Matrice de stockage totale isym=2 IF (IRIGB.EQ.IRIGC) THEN IF (irigd2.irigel(7,irig).eq.0) THEN if (irige2.ne.0) then isym=irige2.irigel(7,irig) else isym=0 endif endif endif call rdscr1(ipt1,des3,isym,ipt2,des2) IF (IERR.NE.0) RETURN irig2.irigel(1,irig)=ipt2 irig2.irigel(3,irig)=des2 IRIG2.IRIGEL(7,irig)=isym * * Remplissage des matrices élémentaires * nelrig=nlrig nligrp=des3.noelep(/1) nligrd=des3.noeled(/1) segini xmatr3 xmatr3.symre=isym xmatrd=irigd2.irigel(4,irig) * segprt,xmatrd xmatrc=irigc2.irigel(4,irig) if (irigb.ne.irigc) then xmatrb=irigb2.irigel(4,irig) else xmatrb=xmatrc endif if (irige.ne.0) then xmatre=irige2.irigel(4,irig) else xmatre=0 endif ntot=nligp+nligd segini amat if (lvdec) then segini amats segini amatv endif segini vvec do ilrig=1,nlrig * Copies do iligp=1,nligp do iligd=1,nligp * segprt,xmatrd * write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd amat(iligd,iligp)=xmatrd.re(iligd,iligp,ilrig) * write(ioimp,*) 'amat(iligd,iligp)=',amat(iligp,iligd) enddo enddo do iligp=1,nligp do iligd=1,nligd * write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd amat(iligd+nligp,iligp)= $ xmatrc.re(iligd,iligp,ilrig) * write(ioimp,*) 'amat(iligd+nligp,iligp)=',amat(iligd * $ +nligp,iligp) enddo enddo do iligp=1,nligp do iligd=1,nligd * write(ioimp,*) 'ilrig,iligp,iligd=',ilrig,iligp,iligd amat(iligp,iligd+nligp)= $ xmatrb.re(iligd,iligp,ilrig) * write(ioimp,*) 'amat(iligp,iligd+nligp)=',amat(iligp * $ ,iligd+nligp) enddo enddo if (xmatre.ne.0) then do iligp=1,nligd do iligd=1,nligd amat(iligp+nligp,iligd+nligp)= $ xmatre.re(iligd,iligp,ilrig) enddo enddo else do iligp=1,nligd do iligd=1,nligd amat(iligp+nligp,iligd+nligp)=0.D0 enddo enddo endif * Sauvegarde de A dans As if (lvdec) then do j=1,ntot do i=1,ntot amats(i,j)=amat(i,j) enddo enddo * segprt,amats endif * if (ilrig.eq.1) segprt,amats * Echelle par la partie haut gauche xscahg=0.d0 do j=1,nligp xscahg=max(xscahg,amat(j,j)) enddo * xpethg positi xpethg=xscahg*xzpiv if (xpethg.EQ.0.d0) then write(ioimp,*) 'D is a null matrix' segprt,des3 write(ioimp,*) 'irig,nrigel=',irig,nrigel write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig write(ioimp,*) 'nligp,nligd=',nligp,nligd segprt,amat return endif * xscabd=0.d0 do j=1,nligp xest=0.d0 xpiv=max(xpethg,amat(j,j)) do i=1,nligd xest=xest-amat(i+nligp,j)*amat(j,i+nligp)/xpiv enddo xscabd=min(xscabd,xest) enddo * xpetbd negatif xpetbd=xscabd*xzpiv if (xpetbd.EQ.0.d0) then write(ioimp,*) 'B or C are null matrices' segprt,des3 write(ioimp,*) 'irig,nrigel=',irig,nrigel write(ioimp,*) 'ilrig,nlrig=',ilrig,nlrig write(ioimp,*) 'nligp,nligd=',nligp,nligd segprt,amat return endif if (lvdec) then xtol=max(xscahg,abs(xscabd))*xzfac endif * * Factorisation de amat (LDMt ou LDLt d'apres Golub et Van Loan 3eme * edition pp. 137-139) * Amélioration possible : prendre en compte la symétrie do j=1,ntot * Résoudre L(1:1j,1:j) v(1:j) = A(1:j,j) do k=1,j vvec(k)=amat(k,j) enddo do k=1,j-1 do i=k+1,j vvec(i)=vvec(i)-vvec(k)*amat(i,k) enddo enddo * Calculer M(j,1:j-1) et stocker dans A (1:j-1,j) do i=1,j-1 amat(i,j)=vvec(i)/amat(i,i) enddo * Stocker d(j) dans A(j,j) amat(j,j)=vvec(j) if (j.le.nligp) then if (abs(amat(j,j)).le.xpethg) then if (ldbg) then write(ioimp,*) 'Pivot petit detecte ilrig,j=' $ ,ilrig,j write(ioimp,*) '| ',amat(j,j),' | < ',xpethg write(ioimp,*) 'Noeud ',ipt1.num(des1.noelep(j) $ ,ilrig), ' Inconnue ',des1.lisinc(j) segprt,amats write(ioimp,*) 'Diag hg :' write(ioimp,*) (amats(jj,jj),jj=1,nligp) write(ioimp,*) 'Diag bd :' write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot) write(ioimp,*) 'Pivots hg :' write(ioimp,*) (amat(jj,jj),jj=1,nligp) write(ioimp,*) 'Pivots bd :' write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot) endif amat(j,j)=xpethg npivhg=npivhg+1 endif else if (abs(amat(j,j)).le.-xpetbd) then if (ldbg) then write(ioimp,*) 'Pivot petit detecte ilrig,j=' $ ,ilrig,j write(ioimp,*) '| ',amat(j,j),' | < -1* ',xpetbd write(ioimp,*) 'Noeud ',ipt1.num(des3.noelep(j $ -nligp),ilrig), ' Inconnue ',des3.lisinc(j $ -nligp) segprt,amats write(ioimp,*) 'Diag hg :' write(ioimp,*) (amats(jj,jj),jj=1,nligp) write(ioimp,*) 'Diag bd :' write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot) write(ioimp,*) 'Pivots hg :' write(ioimp,*) (amat(jj,jj),jj=1,nligp) write(ioimp,*) 'Pivots bd :' write(ioimp,*) (amat(jj,jj),jj=nligp+1,ntot) endif amat(j,j)=xpetbd npivbd=npivbd+1 endif endif * Calculer L(j+1:n,j) et stocker dans A (j+1:n,j) do k=1,j-1 do i=j+1,ntot amat(i,j)=amat(i,j)-vvec(k)*amat(i,k) enddo enddo * write(ioimp,*) 'ilrig,j,a(j,j),xpetbd,xpethg=',ilrig,j * $ ,amat(j,j),xpetbd,xpethg do i=j+1,ntot amat(i,j)=amat(i,j)/amat(j,j) enddo enddo * Vérification de la factorisation if (lvdec) then * segprt,amat do k=1,ntot do i=1,ntot amatv(i,k)=0.D0 enddo enddo do k=1,ntot do i=1,ntot do j=1,min(i,k) xajj=amat(j,j) if (k.eq.j) then xukj=1.d0 else xukj=amat(j,k) endif if (i.eq.j) then xlij=1.d0 else xlij=amat(i,j) endif amatv(i,k)=amatv(i,k)+xlij*xajj*xukj enddo enddo enddo * segprt,amatv do k=1,ntot do i=1,ntot xda=abs(amatv(i,k)-amats(i,k)) if (xda.gt.xtol) then iprint=iprint+1 write(ioimp,*) $ 'Erreur factorisation detectee ilrig=' $ ,ilrig write(ioimp,*) 'aik=',amats(i,k),'aikv=',amatv(i $ ,k) write(ioimp,*) 'i=',i,' k=',k write(ioimp,*) 'xpethg,xpetbd,xtol=',xpethg $ ,xpetbd,xtol write(ioimp,*) 'Diag hg orig :' write(ioimp,*) (amats(jj,jj),jj=1,nligp) write(ioimp,*) 'Diag bd orig :' write(ioimp,*) (amats(jj,jj),jj=nligp+1,ntot) write(ioimp,*) 'Diag hg :' write(ioimp,*) (amatv(jj,jj),jj=1,nligp) write(ioimp,*) 'Diag bd :' write(ioimp,*) (amatv(jj,jj),jj=nligp+1,ntot) segprt,amats segprt,amatv ierr=1 return * write(ioimp,*) 'xda = ',xda,' > ',xtol * write(ioimp,*) 'Noeud dual i ',ipt1.num(des3.noeled(i) * $ ,ilrig), ' Inconnue ',des3.lisdua(i) * write(ioimp,*) 'Noeud primal k',ipt1.num(des3.noelep(k) * $ ,ilrig), ' Inconnue ',des3.lisinc(k) if (iprint.gt.nprint) then return endif endif enddo enddo endif * Produit LDMt avec le bloc en bas à droite et stockage dans RE * Amélioration possible : prendre en compte la symétrie do k=1,nligd k2 = nligp+k do i=1,nligd do j=1,min(i,k) if (k.eq.j) then xukj=1.d0 else endif if (i.eq.j) then xlij=1.d0 else endif xmatr3.re(i,k,ilrig)=xmatr3.re(i,k,ilrig) $ +xlij*xajj*xukj enddo enddo enddo ENDDO segsup vvec if (lvdec) then segsup amats segsup amatv endif segsup amat IRIG2.IRIGEL(4,irig)=xmatr3 IRIG2.IRIGEL(5,irig)=IRIGC2.IRIGEL(5,irig) enddo * Impression pivots nuls npivnu=npivhg+npivbd if (npivnu.gt.0.and.iimpi.ne.0) then * if (npivnu.gt.0) then xpivel=float(npivnu)/float(nelpnu) xpilhg=float(npivhg)/float(nelpnu) xpilbd=float(npivbd)/float(nelpnu) write(ioimp,*) npivnu,' pivots nuls detectes / ',nelpnu, $ ' elements = ',xpivel,' pivots nuls/element' write(ioimp,*) ' dont : ', xpilhg $ ,' pivots nuls bloc(1,1)/element' write(ioimp,*) ' ', xpilbd $ ,' pivots nuls bloc(2,2)/element' endif * segprt,irig2 * * Menage * DO irig=1,nbsou * descr=irigc2.irigel(3,irig) * segsup descr xmatri=irigc2.irigel(4,irig) segsup xmatri IF (IRIGB.NE.IRIGC) THEN * descr=irigb2.irigel(3,irig) * segsup descr xmatri=irigb2.irigel(4,irig) segsup xmatri ENDIF * descr=irigd2.irigel(3,irig) * segsup descr xmatri=irigd2.irigel(4,irig) segsup xmatri IF (IRIGE.NE.0) THEN * descr=irige2.irigel(3,irig) * segsup descr xmatri=irige2.irigel(4,irig) segsup xmatri ENDIF ENDDO segsup irigc2 IF (IRIGB.NE.IRIGC) segsup irigb2 segsup irigd2 IF (IRIGE.NE.0) segsup irige2 * icpr=jcpr inode=jnode jelnum=kelnum izone=jzone segsup izone segsup jelnum segsup inode segsup icpr SEGSUP INFMEL * segsup meleme RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales