c234567 program princip c calcul de la solution de l equation du ptentiel c par elements finis quadratiques Q9 c conditions aux limites : c de Dirichlet sur Gamma_1 u = g c de Neumann sur Gamma_n du/dn = 0 c reservation memoire c234567 implicit real*8(c) implicit real*8(a-b,d-h,o-z) parameter(nmem=1700000) c#include "param" parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0, & h1=1.d0) dimension atab(nmem) dimension itab(nmem) c calcul des dimensions des tableaux c nombre de noeuds nbsom = (nhmai+1)*(nvmai+1) c nombre d'éléments nbel = nhmai/2*nvmai/2 write(*,*)'nombre d elements',nbel write(*,*) 'nombre de noeuds',nbsom c nombre de ddl (1 par noeud) ndl=nbsom c calcul des pointeurs entiers write(*,*) 'calcul des pointeurs' c connectivite elts->ddl (9*nbel) iptie=1 c tableau des noeuds imposes (nbsom) ipticlim=iptie+9*nbel c hauteurs des colonnes (ndl) iptkhj=ipticlim+nbsom c tableau des debuts de colonnes (ndl+1) iptkld=iptkhj+ndl c calcul des pointeurs reels et complexes c (un complexe = 2 reels) c coordonnees des noeuds (2*nbsom) iptx=1 ipty=iptx+nbsom c valeurs des noeuds imposes (ndl) iptclim=ipty+2*nbsom c diagonale de K (ndl) iptvkgd=iptclim+nbsom c triangle sup. de K (lvkgs : inconnu !) iptvkgs=iptvkgd+ndl c lecture du maillage write(*,*) 'lecture du maillage' call lecm(nbel,nbsom,nbclim,itab(iptie), $ atab(iptx),atab(ipty), $ itab(ipticlim),atab(iptclim),xom) c initialisations de : c hauteurs de colonnes c et tables des debuts de colonnes write(*,*) 'initialisations' call init(nbel,nbsom,ndl,itab(iptie) &,itab(iptkhj), & itab(iptkld),lvkgs) c test de memoire write(*,*) 'longueur de vkgs ',lvkgs write(*,*) 'reste',nmem-lvkgs-iptvkgs,' de libre' if (lvkgs+iptvkgs.gt.nmem) then write(*,*) 'pas assez de memoire' write(*,*) 'modifier nmem' stop endif c vecteur second membre (ndl) c234567 iptforce=iptvkgs+lvkgs c assemblage de la matrice de rigidite c assemblage du vecteur second membre c imposition des c-l par penalisation write(*,*) 'assemblage' call assemb(nbclim,nbel,ndl,itab(iptie),itab(iptkld) &, atab(iptvkgd),atab(iptvkgs) &, lvkgs,atab(iptforce),atab(iptx),atab(ipty) &,itab(ipticlim),atab(iptclim),xom) c pointeur de la solution (ndl) iptsolu=iptforce+ndl c factorisation de la matrice de rigidite write(*,*) 'factorisation' call sol(atab(iptvkgs),atab(iptvkgd) &,atab(iptvkgs),atab(iptforce), &itab(iptkld),atab(iptsolu),ndl,6,1,0,0,energ,ier) c calcul de la solution write(*,*) 'resolution' call sol(atab(iptvkgs),atab(iptvkgd) &,atab(iptvkgs),atab(iptforce), & itab(iptkld),atab(iptsolu),ndl,6,0,1,0,energ,ier) c sauvegarde de la solution et calcul de l'erreur c pointeur de la vrai sol (ndl) iptvraisol=iptsolu+ndl write(*,*) 'calcul du gradient' call sauv(nbel,nbsom,itab(iptie), & atab(iptx),atab(ipty),atab(iptsolu) &,atab(iptvraisol)) lastptr=iptvraisol+ndl write(*,*) 'reels : mem total',nmem &,' utilise',lastptr lastpti=iptkld+ndl+1 write(*,*) 'entiers : mem total',nmem &,' utilise',lastpti end c234567 subroutine assemb(nbclim,nel,ndl,ie,kld, & vkgd,vkgs,lvkgs,force,xx,yy, & iclim,clim,xom) c assemblage du second membre et c de la matrice de rigidite implicit real*8(c) implicit real*8(a-b,d-h,o-z) parameter(npg=9,ndle=9) integer ie(ndle,nel),kld(ndl+1) dimension force(ndl),vkgs(lvkgs),vkgd(ndl) &,xx(ndl),yy(ndl) dimension xpg(npg),xref(ndle),yref(ndle) & ,vkelem(ndle,ndle),ypg(npg),hpg(npg), & dmat(3,3),xp(2),vni(ndle,3,npg),ajac(2,2), & ajinv(2,2),bmat(3,ndle),dbmat(3,ndle),btdb(ndle,ndle) &,iclim(ndl),clim(ndl) c element de reference xref(1)=0 xref(2)=.5d0 xref(3)=1.d0 xref(4)=1.d0 xref(5)=1.d0 xref(6)=.5d0 xref(7)=0 xref(8)=0 xref(9)=.5d0 yref(1)=0 yref(2)=0 yref(3)=0 yref(4)=.5d0 yref(5)=1.d0 yref(6)=1.d0 yref(7)=1.d0 yref(8)=.5d0 yref(9)=.5d0 c stockage des points de gauss et des poids xpi=4.d0*datan(1.d0) ab=dsqrt(3.d0/5.d0) xpg(1)=-ab xpg(2)=0.d0 xpg(3)=ab xpg(4)=-ab xpg(5)=0.d0 xpg(6)=ab xpg(7)=-ab xpg(8)=0.d0 xpg(9)=ab ypg(1)=-ab ypg(2)=-ab ypg(3)=-ab ypg(4)=0 ypg(5)=0 ypg(6)=0 ypg(7)=ab ypg(8)=ab ypg(9)=ab do i=1,9 xpg(i)=(xpg(i)+1.d0)*.5d0 ypg(i)=(ypg(i)+1.d0)*.5d0 end do ap=5.d0/9.d0 bp=8.d0/9.d0 cp=5.d0/9.d0 ap=.5d0*ap bp=.5d0*bp cp=.5d0*cp hpg(1)=ap*ap hpg(2)=bp*ap hpg(3)=cp*ap hpg(4)=ap*bp hpg(5)=bp*bp hpg(6)=cp*bp hpg(7)=ap*cp hpg(8)=bp*cp hpg(9)=cp*cp c calcul des fonctions de base c et des derivees aux c points de gauss c234567 un=1.d0 de=2.d0 qa=4.d0 do ipg=1,npg x=xpg(ipg) y=ypg(ipg) c fonction au point de gauss vni(1,1,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)*(y-1.D0) vni(2,1,ipg)= -0.8D1*x*(x-1.D0)*(y-0.5D0)*(y-1.D0) vni(3,1,ipg)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)*(y-1.D0) vni(4,1,ipg)= -0.8D1*x*(x-0.5D0)*y*(y-1.D0) vni(5,1,ipg)= 0.4D1*x*(x-0.5D0)*y*(y-0.5D0) vni(6,1,ipg)= -0.8D1*x*(x-1.D0)*y*(y-0.5D0) vni(7,1,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*y*(y-0.5D0) vni(8,1,ipg)= -0.8D1*(x-0.5D0)*(x-1.D0)*y*(y-1.D0) vni(9,1,ipg)= 0.16D2*x*(x-1.D0)*y*(y-1.D0) c derivee /xi de la fonction au point de gauss vni(1,2,ipg)= 0.4D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(y #-0.5D0)*(y-1.D0) vni(2,2,ipg)= -0.8D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)-0.8D1*x*(y-0.5D0) #*(y-1.D0) vni(3,2,ipg)= 0.4D1*(x-0.5D0)*(y-0.5D0)*(y-1.D0)+0.4D1*x*(y-0.5D0) #*(y-1.D0) vni(4,2,ipg)= -0.8D1*(x-0.5D0)*y*(y-1.D0)-0.8D1*x*y*(y-1.D0) vni(5,2,ipg)= 0.4D1*(x-0.5D0)*y*(y-0.5D0)+0.4D1*x*y*(y-0.5D0) vni(6,2,ipg)= -0.8D1*(x-1.D0)*y*(y-0.5D0)-0.8D1*x*y*(y-0.5D0) vni(7,2,ipg)= 0.4D1*(x-1.D0)*y*(y-0.5D0)+0.4D1*(x-0.5D0) # *y*(y-0.5D0) vni(8,2,ipg)= -0.8D1*(x-1.D0)*y*(y-1.D0)-0.8D1*(x-0.5D0) #*y*(y-1.D0) vni(9,2,ipg)= 0.16D2*(x-1.D0)*y*(y-1.D0)+0.16D2*x*y*(y-1.D0) c derivee /eta de la fonction au point de gauss vni(1,3,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(x #-1.D0)*(y-0.5D0) vni(2,3,ipg)=-0.8D1*x*(x-1.D0)*(y-1.D0)-0.8D1*x*(x-1.D0)*(y-0.5D0) vni(3,3,ipg)=0.4D1*x*(x-0.5D0)*(y-1.D0)+0.4D1*x*(x-0.5D0) #*(y-0.5D0) vni(4,3,ipg)= -0.8D1*x*(x-0.5D0)*(y-1.D0)-0.8D1*x*(x-0.5D0)*y vni(5,3,ipg)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)+0.4D1*x*(x-0.5D0)*y vni(6,3,ipg)= -0.8D1*x*(x-1.D0)*(y-0.5D0)-0.8D1*x*(x-1.D0)*y vni(7,3,ipg)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)+0.4D1*(x-0.5D0)*( #x-1.D0)*y vni(8,3,ipg)= -0.8D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)-0.8D1*(x-0.5D0)*( #x-1.D0)*y vni(9,3,ipg)= 0.16D2*x*(x-1.D0)*(y-1.D0)+0.16D2*x*(x-1.D0)*y enddo c calcul de la matrice D do ii=1,3 do jj=1,3 dmat(ii,jj)=0.d0 enddo enddo dmat(1,1)=xom**2 dmat(2,2)=-1.d0 dmat(3,3)=-1.d0 c boucle sur les elements c assemblage de la matrice de rigidite do ine=1,nel do ii=1,ndle do jj=1,ndle vkelem(ii,jj)=0.d0 enddo enddo c calcul de la matrice de rigidite elementaire do ipg=1,npg c calcul de la matrice jacobienne do ii=1,2 do jj=1,2 ajac(ii,jj)=0.d0 enddo enddo do k=1,ndle xp(1)=xx(ie(k,ine)) xp(2)=yy(ie(k,ine)) c$$$ xp(1)=yref(k)*yref(k) c$$$ xp(2)=-xref(k) do ii=1,2 do jj=1,2 ajac(ii,jj)=ajac(ii,jj)+xp(ii)*vni(k,jj+1,ipg) enddo enddo enddo c$$$ write(*,*) ajac(1,1),' ',ajac(1,2) c$$$ write(*,*) ajac(2,1),' ',ajac(2,2) c$$$ stop c calcul du jacobien det=ajac(1,1)*ajac(2,2)-ajac(1,2)*ajac(2,1) if (det.le.0.d0) then write(*,*) 'det negatif sur l element',ine,det endif c calcul de la jacobienne inverse ajinv(1,1)=ajac(2,2)/det ajinv(2,2)=ajac(1,1)/det ajinv(1,2)=-ajac(1,2)/det ajinv(2,1)=-ajac(2,1)/det c calcul de la matrice B do jj=1,ndle bmat(1,jj)=vni(jj,1,ipg) bmat(2,jj)=0.d0 bmat(3,jj)=0.d0 enddo do j=1,ndle do ii=2,3 do kk=1,2 bmat(ii,j)=bmat(ii,j)+ajinv(kk,ii-1)*vni(j,kk+1,ipg) enddo enddo enddo c calcul de D*B do ii=1,3 do jj=1,ndle dbmat(ii,jj)=0.d0 do kk=1,3 dbmat(ii,jj)=dbmat(ii,jj)+dmat(ii,kk)*bmat(kk,jj) enddo enddo enddo c calcul de Bt*D*B do ii=1,ndle do jj=1,ndle btdb(ii,jj)=0.d0 do kk=1,3 btdb(ii,jj)=btdb(ii,jj)+bmat(kk,ii)*dbmat(kk,jj) enddo enddo enddo do ii=1,ndle do jj=1,ndle vkelem(ii,jj)=vkelem(ii,jj)+btdb(ii,jj)*det*hpg(ipg) enddo enddo c fin de la boucle sur les points de gauss enddo c distribution dans la matrice ligne de ciel do ii=1,ndle do jj=1,ndle i=ie(ii,ine) j=ie(jj,ine) if (i.eq.j) then vkgd(i)=vkgd(i)+vkelem(ii,ii) endif if (i.lt.j) then l=kld(j+1)-j+i vkgs(l)=vkgs(l)+vkelem(ii,jj) endif enddo enddo c fin de la boucle sur les elements enddo c boucle sur les elements c assemblage du second membre do ine=1,nel do ii=1,ndle i=ie(ii,ine) do ipg=1,npg force(i)=0.d0 enddo enddo enddo c imposition des cond. aux limites c par penalisation do ii=1,nbclim i=iclim(ii) vkgd(i)=1.d20 force(i)=clim(ii)*1.d20 enddo return end c234567 subroutine init(nel,nbsom,ndl,ie &,khj, kld,lvkgs) c initialisations de : c hauteurs de colonnes c debuts de colonnes c234567 implicit complex*16(c) implicit real*8(a-b,d-h,o-z) integer ie(9,nel),khj(ndl),kld(ndl+1) c calcul des hauteurs de colonnes do i=1,ndl khj(i)=0 enddo do ine=1,nel do ii=1,9 do jj=1,9 i=ie(ii,ine) j=ie(jj,ine) if(i.lt.j) then if(j-i.gt.khj(j))then khj(j)=j-i endif endif enddo enddo enddo do i=1,ndl c write(*,*) 'hauteur',i,' =',khj(i) enddo c calcul de kld kld(1)=1 do i=2,ndl+1 kld(i)=kld(i-1)+khj(i-1) enddo c calcul dela longueur de vkgs lvkgs=kld(ndl+1)-1 do i=1,ndl+1 c write(*,*) 'debut col.',i,' =',kld(i) enddo return end subroutine lecm(nbel,nbsom,nbclim,ie,xx,yy,iclim,clim,xom) c lecture du maillage et de la connectivite implicit real*8(c) implicit real*8(a-b,d-h,o-z) c dimension : fichier tanaka, nbre de points horizontaux et verticaux c#include "param" parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0, & h1=1.d0) dimension xx(nbsom),yy(nbsom),clim(nbsom) integer kinterm(16) dimension ie(9,nbel),iclim(nbsom) dimension xtanak(ntanak),ytanak(ntanak),faitanak(ntanak) dimension xmaiup(0:nhmai),ymaiup(0:nhmai),faiup(0:nhmai) c ouverture du fichier des valeurs données par c la méthode de tanaka open(1,file='tanaka.dat',status='old') do i=1,ntanak read(1,*) j,xtanak(i),ytanak(i),faitanak(i) end do close(1) c nombre de noeuds nbsom = (nhmai+1)*(nvmai+1) c nombre d'éléments nbel = nhmai/2*nvmai/2 c write(*,*) 'nbre de sommets=',nbsom c write(*,*) 'nbre d elements=',nbel c longueur du domaine de calcul (adim) c along = 13.d0 c hauteur d'eau à l'infini (adim) c h1=1.d0 c pas du maillage en x dx = along/nhmai do i = 0,nhmai xmaiup(i)=i*dx c recherche de l'intervalle ou se trouve le point j=0 555 j=j+1 if (xmaiup(i).ge.xtanak(j)) then goto 555 endif c interpolation linéaire entre j-1 et j t=(xmaiup(i)-xtanak(j-1))/(xtanak(j)-xtanak(j-1)) ymaiup(i)=t*ytanak(j)+(1-t)*ytanak(j-1) faiup(i)=t*faitanak(j)+(1-t)*faitanak(j-1) end do c$$$ do i = 0,nhmai c$$$ write(*,*) i,xmaiup(i),ymaiup(i),faiup(i) c$$$ enddo c$$$ stop c calcul des coor des sommets do i =0,nhmai do j =0,nvmai dy=(ymaiup(i)+h1)/nvmai ii=j+(nvmai+1)*i+1 xx(ii)=xmaiup(i) yy(ii)=-h1+j*dy end do end do c$$$ open(2,file='g1',status='unknown') c$$$ do i=1,nbsom c$$$ write(2,*) xx(i),yy(i) c$$$ enddo c$$$ close(2) c calcul de la connectivite do i=0,nhmai-1,2 do j=0,nvmai-1,2 iel = j/2+(nvmai/2)*i/2+1 ii1=(j+0)+(nvmai+1)*(i+0)+1 ii2=(j+0)+(nvmai+1)*(i+1)+1 ii3=(j+0)+(nvmai+1)*(i+2)+1 ii4=(j+1)+(nvmai+1)*(i+2)+1 ii5=(j+2)+(nvmai+1)*(i+2)+1 ii6=(j+2)+(nvmai+1)*(i+1)+1 ii7=(j+2)+(nvmai+1)*(i+0)+1 ii8=(j+1)+(nvmai+1)*(i+0)+1 ii9=(j+1)+(nvmai+1)*(i+1)+1 ie(1,iel)=ii1 ie(2,iel)=ii2 ie(3,iel)=ii3 ie(4,iel)=ii4 ie(5,iel)=ii5 ie(6,iel)=ii6 ie(7,iel)=ii7 ie(8,iel)=ii8 ie(9,iel)=ii9 end do end do c$$$ open(3,file='g2',status='unknown') c$$$ do i=1,nbel c$$$ write(3,*) xx(ie(1,i)),yy(ie(1,i)) c$$$ write(3,*) xx(ie(3,i)),yy(ie(3,i)) c$$$ write(3,*) xx(ie(5,i)),yy(ie(5,i)) c$$$ write(3,*) xx(ie(7,i)),yy(ie(7,i)) c$$$ enddo c$$$ close(3) c$$$ stop c$$$ do i=1,nbel c$$$ write(*,*) i,'-->',ie(1,i),ie(2,i),ie(3,i),ie(4,i) c$$$ & ,ie(5,i),ie(6,i),ie(7,i),ie(8,i),ie(9,i) c$$$ enddo c$$$ stop c calcul des conditions aux limites c bord gauche nbclim=0 do j=0,nvmai-1 nbclim=nbclim+1 i=0 ii=j+(nvmai+1)*i+1 clim(nbclim)=0 c clim(nbclim)=1 iclim(nbclim)=ii end do c partie sup do i=0,nhmai nbclim=nbclim+1 j=nvmai ii=j+(nvmai+1)*i+1 clim(nbclim)=faiup(i) c clim(nbclim)=1 iclim(nbclim)=ii end do c affichage xom=0 write(*,*) 'nbre de c_l=',nbclim c$$$ do i=1,nbclim c$$$ xa=xx(iclim(i)) c$$$ ya=yy(iclim(i)) c$$$ write(*,*) i,' no sommet=',iclim(i),' val=',clim(i) c$$$ enddo return end subroutine sauv(nbel,nbsom,ie,xx,yy,sol,vraisol) c lecture du maillage et de la connectivite implicit real*8(c) implicit real*8(a-b,d-h,o-z) c#include "param" parameter(ntanak=80,nhmai=260,nvmai=60,along=13.d0, & h1=1d0) dimension xx(nbsom),yy(nbsom) dimension ie(9,nbel),sol(nbsom), & vraisol(nbsom),vni(9,3),ajac(2,2), & ajinv(2,2),vnivrai(9,3) dimension xref(9),yref(9) c ouverture des fichiers c$$$ open(1,file='MOSMEF.SOL',status='unknown') c$$$ c$$$ write(1,*) c$$$ write(1,*) c$$$ write(1,'(6i5)') 1,1,1,0,0,nbsom c$$$ write(1,*)' CAS NO : 0 PAS NUMERO : 0 ', c$$$ & 'ITER. NO : 0 0.00000E+00' c$$$ c$$$c ecriture de la solution c$$$ c$$$ xpi=4.d0*datan(1.d0) c$$$ xk=dsqrt(2.d0)/2.d0 c$$$ yk=dsqrt(2.d0)/2.d0 c$$$ c$$$ do is=1,nbsom c$$$ c$$$c write(*,*) is,xx(is),yy(is) c$$$ xa=xx(is) c$$$ ya=yy(is) c$$$c val=sol(is)-dsin(2.d0*xpi*(xk*xa+yk*ya)) c$$$ val=sol(is) c$$$ c$$$ write(1,'(2i5,4e12.5)') is,4,xx(is),yy(is),0.,val c$$$ c$$$ enddo c$$$ c$$$ close(1) c calcul du gradient du potentiel aux points du fichier c 'pcal.dat' c longueur du domaine de calcul (adim) c along = 13.d0 c hauteur d'eau à l'infini (adim) c h1=1.d0 c pas du maillage en x pasx = along/nhmai c distance et vitesse caracteristiques dchar=0.310d0 vchar=dsqrt(9.81d0*dchar) open(17,file='pcal.dat',status='old') open(18,file='plotuv.dat',status='unknown') c nombre de points où calculer le gradient read(17,*) nuser c write(18,*) nuser do iu=1,nuser read(17,*) if,xuser,yuser c adimensionnement xuser=xuser/dchar yuser=yuser/dchar c write(*,*) xuser,yuser c recadrage en abscisse xsauv=xuser ysauv=yuser if (xuser.lt.0) then xuser = -xuser c write(*,*) xuser,iu c stop endif if (xuser.gt.along) then xuser = along-1.d-8 endif if (yuser.lt.-1) then write(*,*) 'erreur pt n°',iu,': y=',yuser*dchar,dchar stop endif c recherche de l'élément dans lequel se trouve le point i=int(xuser/2.d0/pasx) kuser = 0 do j=0,nvmai/2-1 k= j+(nvmai/2)*i+1 if (k.gt.nbel) stop x1=xx(ie(1,k)) y1=yy(ie(1,k)) x2=xx(ie(3,k)) y2=yy(ie(3,k)) x3=xx(ie(5,k)) y3=yy(ie(5,k)) x4=xx(ie(7,k)) y4=yy(ie(7,k)) t=(xuser-x1)/(x2-x1) yb= t*y2+(1-t)*y1 yh= t*y3+(1-t)*y4 c write(*,*) k,x1,y1,x2,y2,x3,y3,x4,y4 if (yuser.le.yh.and.yuser.ge.yb) then kuser = k c write(*,*) k endif end do if (kuser.eq.0) then write(*,*) 'pt n°',iu,' dans l''air ', & 'x=',xsauv*dchar,' y=',ysauv*dchar endif if (kuser.ne.0) then write(*,*) 'pt n°',iu,' dans l''eau ', & 'x=',xsauv*dchar,' y=',ysauv*dchar hmax=yh c write(*,*) hmax c recherche des coordonnees dans l elt de ref c par newton x = .5d0 y = .5d0 c image par tau de xchap,ychap c fonction de base 111 continue c write(*,*) x,y,det,kuser call fbase(vni,x,y,ajac,ajinv,det,vnivrai,ie,xx,yy,kuser) taux=0 tauy=0 do in=1,9 taux=taux+vni(in,1)*xx(ie(in,kuser)) tauy=tauy+vni(in,1)*yy(ie(in,kuser)) end do dx = ajinv(1,1)*(xuser-taux)+ajinv(1,2)*(yuser-tauy) dy = ajinv(2,1)*(xuser-taux)+ajinv(2,2)*(yuser-tauy) x=x+dx y=y+dy if (dsqrt(dx*dx+dy*dy).ge.1.d-10) goto 111 c test if (x.gt.1.01d0.or.y.gt.1.01d0 $ .or.x.lt.-0.01d0.or.y.lt.-0.01d0) then write(*,*) x,y stop endif c calcul de phi valphi=0; gx=0; gy=0; do in=1,9 valphi=valphi+vni(in,1)*sol(ie(in,kuser)) gx=gx+vnivrai(in,2)*sol(ie(in,kuser)) gy=gy+vnivrai(in,3)*sol(ie(in,kuser)) end do else valphi=0 gx=0 gy=0 endif if (xsauv.lt.0) then valphi=-valphi gy=-gy endif c calcul de la pression par Bernouilli xpinf=1.d5 rho=1000 gpes=9.81 c vref=dsqrt(gpes*h1) vref=1 hinf=yy(ie(7,nbel)) xpress = xpinf +rho*gpes*(hmax-yuser) c $ -rho*0.5d0*(gx*gx+gy*gy)*vref**2 xpress = xpress if (kuser.eq.0) then xpress = 1d5 rho =1 endif write(18,*) iu,xsauv*dchar,ysauv*dchar $ ,rho,gx*vchar,gy*vchar,xpress end do close(17) close(18) c tracés tecplot pour vérifier c$$$ open(19,file='soliton.dat',status='unknown') c$$$ c$$$c element de reference c$$$ c$$$ xref(1)=0 c$$$ xref(2)=.5d0 c$$$ xref(3)=1.d0 c$$$ xref(4)=1.d0 c$$$ xref(5)=1.d0 c$$$ xref(6)=.5d0 c$$$ xref(7)=0 c$$$ xref(8)=0 c$$$ xref(9)=.5d0 c$$$ c$$$ yref(1)=0 c$$$ yref(2)=0 c$$$ yref(3)=0 c$$$ yref(4)=.5d0 c$$$ yref(5)=1.d0 c$$$ yref(6)=1.d0 c$$$ yref(7)=1.d0 c$$$ yref(8)=.5d0 c$$$ yref(9)=.5d0 c$$$ c$$$ do iel=1,nbel c$$$ do ii=1,7,2 c$$$ c$$$ call fbase(vni,xref(ii),yref(ii),ajac, c$$$ $ ajinv,det,vnivrai,ie,xx,yy,iel) c$$$ gx=0 c$$$ gy=0 c$$$ do in=1,9 c$$$ gx=gx+vnivrai(in,2)*sol(ie(in,iel)) c$$$ gy=gy+vnivrai(in,3)*sol(ie(in,iel)) c$$$ end do c$$$ c$$$ c$$$ write(19,*) xx(ie(ii,iel)),yy(ie(ii,iel)), c$$$ $ sol(ie(ii,iel)),gx,gy c$$$ enddo c$$$ enddo c$$$ c$$$ do iel=1,nbel c$$$ write(19,*) (iel-1)*4+1,(iel-1)*4+2,(iel-1)*4+3,(iel-1)*4+4 c$$$ enddo c$$$ c$$$ close(19) return end subroutine sol (vkgs,vkgd,vkgi,vfg,kld,vu,neq,mp,ifac,isol, $ nsym,energ,ier) c resolution d'un systeme lineaire symetrique ou non. la matrice est c stockee par ligne de ciel,en memoire dans les tables vkgs,vkgd,vkgi c c entrees c vkgs,vkgd,vkgi matrice du systeme : parties superieure, c diagonale, inferieure (double precision) c vfg second membre c kld pointeurs vers les hauts de colonne c vu vecteur solution (qui peut etre vfg) c neq nombre d'equations c mp unite logique d'impression c ifac si ifac.eq.1 triangularisation de c la matrice c isol si isol.eq.1 calcul de la solution a c partir de la matrice triangularisee c nsym indice de probleme non symetrique c sorties c vkgs,vkgd,vkgi matrice triangularisee (si ifac.eq.1) c vfg solution (si isol.eq.1) c energ energie du systeme (si nsym.eq.0) c ier mis a 1 si pivot nul rencontre c c=========================== debut des declarations ==================== implicit real*8 (a-h,o,q-z) implicit integer (p) dimension vkgs(1),vkgd(1),vkgi(1),vfg(1),kld(1),vu(1) data vzero/0.d0/ c=========================== debut du code executable ================== c c------- traitement c ik=1 if(vkgd(1).eq.vzero) goto 800 energ=vzero ier=0 if (isol.eq.1) then do i = 1, neq vu(i) = vfg(i) end do end if c c------- pour chaque colonne ik a modifier c jhk=1 do 100 ik=2,neq c c------- pointeur du haut de la colonne suivante ik+1 c jhk1=kld(ik+1) c c------- hauteur de la colonne ik (hors termes superieur et diagonal) c lhk=jhk1-jhk lhk1=lhk-1 c c------- ligne du premier terme a modifier dans la colonne ik c imin=ik-lhk1 imin1=imin-1 c c------- ligne du dernier terme a modifier dans la colonne ik c imax=ik-1 if(lhk1.lt.0) goto 100 if(ifac.ne.1) goto 90 if(nsym.eq.1) vkgi(jhk)=vkgi(jhk)/vkgd(imin1) if(lhk1.eq.0) goto 40 c c------- modifier les termes non diagonaux de la colonne ik c jck=jhk+1 jhj=kld(imin) c c------- pour chaque terme place en jck, correspondant a la colonne ij c do 30 ij=imin,imax jhj1=kld(ij+1) c c------- nombre de termes modificatifs du terme place en jck c ic=min0(jck-jhk,jhj1-jhj) if(ic.le.0.and.nsym.eq.0) goto 20 c1=vzero if(ic.le.0) goto 17 j1=jhj1-ic j2=jck-ic if(nsym.eq.1) goto 15 vkgs(jck)=vkgs(jck)-scal(vkgs(j1),vkgs(j2),ic) goto 20 15 vkgs(jck)=vkgs(jck)-scal(vkgi(j1),vkgs(j2),ic) c1=scal(vkgs(j1),vkgi(j2),ic) 17 vkgi(jck)=(vkgi(jck)-c1)/vkgd(ij) 20 jck=jck+1 30 jhj=jhj1 c c------- modifier le terme diagonal c 40 jck=jhk cdiag=vzero do 70 ij=imin1,imax c1=vkgs(jck) if(nsym.eq.1) goto 50 c2=c1/vkgd(ij) vkgs(jck)=c2 goto 60 50 c2=vkgi(jck) 60 cdiag=cdiag+c1*c2 70 jck=jck+1 vkgd(ik)=vkgd(ik)-cdiag if (vkgd(ik).eq.0.) goto 800 c c------- resolution du systeme triangulaire inferieur c 90 if(isol.ne.1) goto 100 if(nsym.ne.1) vu(ik)=vfg(ik)-scal(vkgs(jhk),vu(imin1),lhk) if(nsym.eq.1) vu(ik)=vfg(ik)-scal(vkgi(jhk),vu(imin1),lhk) 100 jhk=jhk1 if(isol.ne.1) goto 9999 c c------- resolution du systeme diagonal : c if(nsym.eq.1) goto 120 do 110 ik=1,neq c1=vkgd(ik) if (c1.eq.vzero) goto 800 c2=vu(ik)/c1 vu(ik)=c2 110 energ=energ+c1*c2*c2 c c------- resolution du systeme triangulaire superieur c 120 ik=neq+1 jhk1=kld(ik) 130 ik=ik-1 if(nsym.eq.1) vu(ik)=vu(ik)/vkgd(ik) if(ik.eq.1) goto 9999 c1=vu(ik) jhk=kld(ik) jbk=jhk1-1 if(jhk.gt.jbk)goto 150 ij=ik-jbk+jhk-1 do 140 jck=jhk,jbk vu(ij)=vu(ij)-vkgs(jck)*c1 140 ij=ij+1 150 jhk1=jhk goto 130 c c------- erreurs c 800 if (mp.ne.0) write(mp,8000) ik 8000 format(' * sol pivot nul equation',i5) ier=1 goto 9999 c c------- fin c 9999 continue return c=========================== fin du module sol ================== end function scal(x,y,n) c======================================================================= c calcul du produit scalaire c======================================================================= implicit real*8 (a-h,o-z) dimension x(1),y(1) data zero/0.d0/ c----------------------------------------------------------------------- scal=zero do i=1,n scal=scal+x(i)*y(i) enddo return end subroutine fbase(vni,x,y,ajac,ajinv,det,vnivrai,ie,xx,yy,k) implicit real*8(c) implicit real*8(a-b,d-h,o-z) dimension vni(9,3),ajac(2,2),ajinv(2,2) dimension vnivrai(9,3),ie(9,*),xx(*),yy(*) vni(1,1)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)*(y-1.D0) vni(2,1)= -0.8D1*x*(x-1.D0)*(y-0.5D0)*(y-1.D0) vni(3,1)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)*(y-1.D0) vni(4,1)= -0.8D1*x*(x-0.5D0)*y*(y-1.D0) vni(5,1)= 0.4D1*x*(x-0.5D0)*y*(y-0.5D0) vni(6,1)= -0.8D1*x*(x-1.D0)*y*(y-0.5D0) vni(7,1)= 0.4D1*(x-0.5D0)*(x-1.D0)*y*(y-0.5D0) vni(8,1)= -0.8D1*(x-0.5D0)*(x-1.D0)*y*(y-1.D0) vni(9,1)= 0.16D2*x*(x-1.D0)*y*(y-1.D0) c derivee /xi de la fonction au point de gauss vni(1,2)= 0.4D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(y #-0.5D0)*(y-1.D0) vni(2,2)= -0.8D1*(x-1.D0)*(y-0.5D0)*(y-1.D0)-0.8D1*x*(y-0.5D0) #*(y-1.D0) vni(3,2)= 0.4D1*(x-0.5D0)*(y-0.5D0)*(y-1.D0)+0.4D1*x*(y-0.5D0) #*(y-1.D0) vni(4,2)= -0.8D1*(x-0.5D0)*y*(y-1.D0)-0.8D1*x*y*(y-1.D0) vni(5,2)= 0.4D1*(x-0.5D0)*y*(y-0.5D0)+0.4D1*x*y*(y-0.5D0) vni(6,2)= -0.8D1*(x-1.D0)*y*(y-0.5D0)-0.8D1*x*y*(y-0.5D0) vni(7,2)= 0.4D1*(x-1.D0)*y*(y-0.5D0)+0.4D1*(x-0.5D0) # *y*(y-0.5D0) vni(8,2)= -0.8D1*(x-1.D0)*y*(y-1.D0)-0.8D1*(x-0.5D0) #*y*(y-1.D0) vni(9,2)= 0.16D2*(x-1.D0)*y*(y-1.D0)+0.16D2*x*y*(y-1.D0) c derivee /eta de la fonction au point de gauss vni(1,3)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)+0.4D1*(x-0.5D0)*(x #-1.D0)*(y-0.5D0) vni(2,3)=-0.8D1*x*(x-1.D0)*(y-1.D0)-0.8D1*x*(x-1.D0)*(y-0.5D0) vni(3,3)=0.4D1*x*(x-0.5D0)*(y-1.D0)+0.4D1*x*(x-0.5D0) #*(y-0.5D0) vni(4,3)= -0.8D1*x*(x-0.5D0)*(y-1.D0)-0.8D1*x*(x-0.5D0)*y vni(5,3)= 0.4D1*x*(x-0.5D0)*(y-0.5D0)+0.4D1*x*(x-0.5D0)*y vni(6,3)= -0.8D1*x*(x-1.D0)*(y-0.5D0)-0.8D1*x*(x-1.D0)*y vni(7,3)= 0.4D1*(x-0.5D0)*(x-1.D0)*(y-0.5D0)+0.4D1*(x-0.5D0)*( #x-1.D0)*y vni(8,3)= -0.8D1*(x-0.5D0)*(x-1.D0)*(y-1.D0)-0.8D1*(x-0.5D0)*( #x-1.D0)*y vni(9,3)= 0.16D2*x*(x-1.D0)*(y-1.D0)+0.16D2*x*(x-1.D0)*y c calcul de la jacobienne ajac(1,1)=0 ajac(1,2)=0 ajac(2,1)=0 ajac(2,2)=0 do in=1,9 ajac(1,1)=ajac(1,1)+vni(in,2)*xx(ie(in,k)) ajac(1,2)=ajac(1,2)+vni(in,3)*xx(ie(in,k)) ajac(2,1)=ajac(2,1)+vni(in,2)*yy(ie(in,k)) ajac(2,2)=ajac(2,2)+vni(in,3)*yy(ie(in,k)) end do det=ajac(1,1)*ajac(2,2)-ajac(1,2)*ajac(2,1) if (det.le.0.d0) then write(*,*) 'det negatif',det stop endif c calcul de la jacobienne inverse ajinv(1,1)=ajac(2,2)/det ajinv(2,2)=ajac(1,1)/det ajinv(1,2)=-ajac(1,2)/det ajinv(2,1)=-ajac(2,1)/det c calcul du gradient dans les bonnes variables do ii=1,9 vnivrai(ii,1)=vni(ii,1) vnivrai(ii,2)= ajinv(1,1) * vni(ii,2) + ajinv(2,1) * vni(ii,3) vnivrai(ii,3)= ajinv(1,2) * vni(ii,2) + ajinv(2,2) * vni(ii,3) enddo return end