C BARYC5 SOURCE GOUNAND 26/06/09 21:15:01 12566 SUBROUTINE BARYC5(MELEMX,NKPVIR,TRAVK,KGRAV) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER(I-N) C*********************************************************************** C NOM : BARYC5 C DESCRIPTION : C CALCUL LE BARYCENTRE D'UN MAILLAGE ET LA METRIQUE MOYENNE ASSOCIEE C Repris de baryce.eso. C Par rapport à baryc2, ignore le noeud virtuel KPVIRT C Gère le noeud virtuel KPVIRT C C Par rapport à baryc3, on gère le MCOORD différemment, en vue de C faire moins de SEGADJ à l'aide du segment TRAVK C C Comme baryc4 mais avec un MELEMX C C KGRAV est le numéro du nouveau noeud créé C C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : C APPELES (E/S) : C APPELES (BLAS) : C APPELES (CALCUL) : C APPELE PAR : C*********************************************************************** C SYNTAXE GIBIANE : C ENTREES : C ENTREES/SORTIES : C SORTIES : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 26/09/2017, version initiale C HISTORIQUE : v1, 26/09/2017, création C HISTORIQUE : C HISTORIQUE : C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC CCREEL -INC TMATOP2 -INC TMATOP1 *-INC SMELEMX -INC SMCOORD POINTEUR KCOORD.MCOORD *-INC SMETRIQ POINTEUR KCMETR.METRIQ *-INC STRAVJ POINTEUR TRAVK.TRAVJ LOGICAL lchang PARAMETER(NGRAV=3) DIMENSION XGRAV(NGRAV) PARAMETER(NMET=6) DIMENSION XMET(NMET) DIMENSION A(3,3) DIMENSION IDXSYM(3,3,3) DIMENSION INVSYM(2,6,3) DIMENSION XSUM(3),TLNSUM(3,3) DIMENSION DXGRAV(3) DIMENSION DLMGRA(3,3) LOGICAL LDBG * DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/ DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/ DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/ DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/ * DATA (((INVSYM(I,J,K),I=1,2),J=1,1),K=1,1) /1,1/ DATA (((INVSYM(I,J,K),I=1,2),J=1,3),K=2,2) /1,1,2,1,2,2/ DATA (((INVSYM(I,J,K),I=1,2),J=1,6),K=3,3) $ /1,1,2,1,2,2,3,1,3,2,3,3/ * * * Executable statements * * SEGACT,MCOORD * LDBG=.TRUE. * LDBG=impr.ge.3 LDBG=.FALSE. IDIMP1=IDIM+1 KCOORD=TRAVK.COORD KCMETR=TRAVK.CMETR * mis dans topv2 SEGACT KCOORD*MOD DO I=1,NGRAV XGRAV(I)=0.D0 ENDDO DO I=1,NMET XMET(I)=0.D0 ENDDO NPOIN=0 * SEGACT,MELEMX * DO 3 J=1,IPT1.NUM(/2) * DO 31 I=1,IPT1.NUM(/1) IF (impr.ge.4) write(ioimp,'(A,4I3)') $ 'baryc5 : IMET,NLCOU,NNCOU,NKPVIR=',IMET,NLCOU,NNCOU,NKPVIR DO 3 ILCOU=1,NLCOU DO 31 INCOU=1,NNCOU INO=MELEMX.NUMX(INCOU,ILCOU) IF (NKPVIR.GT.0) THEN * ldbg=.true. IF (INO.LE.NKPVIR) GOTO 31 ENDIF NPOIN=NPOIN+1 IREF=IDIMP1*(INO-1) DO 5 L=1,IDIM XGRAV(L)=XGRAV(L)+KCOORD.XCOOR(IREF+L) 5 CONTINUE IF (KCMETR.NE.0) THEN if (ldbg) then write(ioimp,*) 'ilcou,incou,ino=',ilcou,incou,ino write(ioimp,*) ' coo=',(KCOORD.XCOOR(IREF+L),L=1 $ ,IDIM) write(ioimp,*) ' met=',(KCMETR.XIN(L,INO),L=1 $ ,KCMETR.XIN(/1)) endif DO 6 ININ=1,KCMETR.XIN(/1) XMET(ININ)=XMET(ININ)+KCMETR.XIN(ININ,INO) 6 CONTINUE ENDIF 31 CONTINUE 3 CONTINUE if (ldbg) write(ioimp,*) 'NPOIN=',NPOIN * SEGDES,MELEMX C ON MET LE CENTRE DE GRAVITE DANS LA TABLE DES POINTS C ET LA METRIQUE ASSOCIEE LE CAS ECHANT NPCOUN=TRAVK.NPCOU+1 CALL TOPADP(TRAVK,NPCOUN,1,lchang,'baryc5 : TRAVK') if (ierr.ne.0) return if (iveri.ge.2) then call vetopi(travk,'baryc5 : Apres extension travk 1') if (ierr.ne.0) return endif do i=1,idim xgrav(i)=xgrav(i)/NPOIN enddo IF (KCMETR.NE.0) THEN DO 12 ININ=1,KCMETR.XIN(/1) KCMETR.XIN(ININ,NPCOUN)=XMET(ININ)/NPOIN 12 CONTINUE ENDIF if (ldbg) then write(ioimp,*) 'Centre de gravite' write(ioimp,*) ' coo=',(xgrav(L),L=1,IDIM) write(ioimp,*) ' met=',(KCMETR.XIN(L,NPCOUN),L=1,KCMETR.XIN( $ /1)) endif * * On essaie de mieux placer le point milieu * if (kcmetr.ne.0.and.imomet.eq.1.and.imobar.eq.1) then * tsum=0.d0 do j=1,idim xsum(j)=0.d0 dxgrav(j)=0.d0 do i=1,idim dlmgra(i,j)=0.d0 tlnsum(i,j)=0.d0 enddo enddo if (imet.eq.4) then imx=idim else imx=1 endif DO 4 INL=1,NLCOU DO 41 INN=1,NNCOU INO=MELEMX.NUMX(INN,INL) IF (NKPVIR.GT.0) THEN IF (INO.LE.NKPVIR) GOTO 41 ENDIF IREF=IDIMP1*(INO-1) * +1 car c'est le log de la metrique inverse DO J=1,imx DO I=1,imx * A(I,J)=KCMETR.XIN(IDXSYM(I,J,imx),INO) A(I,J)=(KCMETR.XIN(IDXSYM(I,J,imx),INO) $ +KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN))/2 * A(I,J)=KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN) enddo ENDDO if (ldbg) then write(ioimp,*) 'Log Tenseur metrique inverse voulu' write(ioimp,*) ' a=',((a(L,K),L=1,imx),K=1,imx) endif * Projection suivant la direction centre de gravite-noeud puis exponentielle if (imet.eq.4) then aproj=0.d0 do j=1,idim do i=1,idim aproj=aproj+(KCOORD.XCOOR(IREF+i)-XGRAV(i))*a(i $ ,j)*(KCOORD.XCOOR(IREF+j)-XGRAV(j)) enddo enddo dx2=0.d0 do i=1,idim dx2=dx2+(KCOORD.XCOOR(IREF+i)-XGRAV(i))**2 enddo aproj=aproj/dx2 xmii=exp(aproj) * Imet = 3 ELSE xmii=exp(a(1,1)) ENDIF if (ldbg) then write(ioimp,*) 'Tenseur metrique inverse voulu' write(ioimp,*) ' xmii=',xmii endif tsum=tsum+xmii DO J=1,imx DO I=1,imx tlnsum(i,j)=tlnsum(i,j)+xmii*(KCMETR.XIN(IDXSYM(I,J $ ,imx),INO)-KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN)) ENDDO ENDDO do i=1,idim xsum(i)=xsum(i)+xmii*(KCOORD.XCOOR(IREF+i) $ -XGRAV(i)) enddo 41 CONTINUE 4 CONTINUE DO J=1,imx DO I=1,imx tlnsum(i,j)=tlnsum(i,j)/NPOIN ENDDO ENDDO do i=1,idim xsum(i)=xsum(i)/NPOIN enddo tsum=tsum/NPOIN if (ldbg) then write(ioimp,*) 'Moyennes' write(ioimp,*) ' xsum=',(xsum(L),L=1,IDIM) write(ioimp,*) ' tsum=',tsum write(ioimp,*) ' tlnsum=',((tlnsum(L,K),L=1,imx),K=1,imx) endif xrii=1.d0/tsum DO J=1,imx DO I=1,imx dlmgra(i,j)=dlmgra(i,j)+xrii*tlnsum(i,j) ENDDO ENDDO do i=1,idim dxgrav(i)=dxgrav(i)+xrii*xsum(i) enddo if (ldbg) then write(ioimp,*) 'Depl gravite ameliore ?' write(ioimp,*) ' dxgrav=',(dxgrav(L),L=1,IDIM) write(ioimp,*) 'DLog metrique gravite ameliore ?' write(ioimp,*) ' dlmgra=',((dlmgra(L,K),L=1,imx),K=1,imx) endif do l=1,idim xgrav(l)=xgrav(l)+dxgrav(l) enddo DO ININ=1,KCMETR.XIN(/1) I=INVSYM(1,ININ,IDIM) J=INVSYM(2,ININ,IDIM) KCMETR.XIN(ININ,NPCOUN)=KCMETR.XIN(ININ,NPCOUN)+DLMGRA(I,J) ENDDO * if (ldbg) then write(ioimp,*) 'Centre de gravite ameliore ?' write(ioimp,*) ' coo=',(xgrav(L),L=1,IDIM) write(ioimp,*) ' met=',(KCMETR.XIN(L,NPCOUN),L=1 $ ,KCMETR.XIN(/1)) endif endif * write(ioimp,*) 'npcoun,npcou,npmax',npcoun,travk.npcou,travk.npmax * NBPTS=XCOOR(/1)/IDIMP1+1 * SEGADJ,MCOORD * IREF=(NBPTS-1)*IDIMP1 IREF=(NPCOUN-1)*IDIMP1 DO 11 I=1,IDIMP1 KCOORD.XCOOR(IREF+I)=XGRAV(I) 11 CONTINUE * KGRAV=XCOOR(/1)/IDIMP1 * if (ldbg) then * write(ioimp,*) 'Stop !' * ierr=1 * return * endif KGRAV=NPCOUN if (iveri.ge.2.and.lchang) then call vetopi(travk,'baryc5 : Apres extension travk 2') if (ierr.ne.0) return endif RETURN END