subroutine poiss(pri,isize,ierr) implicit double precision (a-h,o-z) dimension pri(*) character*3 refchr * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact common /dbasco/ ifree,ilast common /pointer/ larhd,lacond,laconv,laabso,lasegd,laisgd * ********************************************************* * initiate lifo-stack ********************************************************* * call dimpr(pri,isize,idum,idum,idum,2) isla=0 ierr=0 imesh=0 * ********************************************************* * set variables ********************************************************* * infile=8 maxt=isize/45 maxv=isize/30 * open(unit=infile,file='area.scratch',status='old') open(unit=77,file='read.scratch',access='sequential', & form='unformatted',status='unknown') open(unit=7,file='plotsol',status='unknown') open(unit=78,file='solu.scratch',status='unknown') open(unit=79,file='solu2.scratch',status='unknown') * * Process input data * call dimpr(pri,3,maxt,4,lanode,1) call dimpr(pri,2,maxv,8,lacoeg,1) call dimpr(pri,maxv,1,4,lanoco,1) call dimpr(pri,3,maxt,4,lanelc,1) call dimpr(pri,4,maxt,8,lavnor,1) call dimpr(pri,3,maxt,4,laited,1) * call inputpoi(pri,infile,6,ierr,pri(lacoeg),pri(lanode), & pri(lanoco),pri(lanelc),refchr) * call dwmatg(pri(lanelc),3,nrtri,0,'i','nelcon') **************** * TEMPORARY **************** * call fooread(79,pri(lavnor),pri(lanelc),nvface) isw=1 * ********************************************************* * Create logical vectors ********************************************************* * call dimpr(pri,maxt,1,4,laleve,1) call dimpr(pri,5,maxt,4,lanelg,1) call dimpr(pri,maxt,1,4,lairef,1) call dimpr(pri,maxv,1,8,lau,1) call dimpr(pri,4,maxt,4,laitri,1) * ********************************************************* * Generate logical arrays ********************************************************* * call gener(pri(lanelg),pri(lanoco),pri(laleve)) * tolref=tol/100.d0 tol=1.d0 * ********************************************************* * tolref is the tolerance we are aiming for * tol is the tolerance we are using to (try to) * achieve equidistributed error * on the current level ********************************************************* * nvact=nno * ********************************************************* * Main loop for refined meshes ********************************************************* * ||||||||||||||||||||||||||||||||||||||||||||||||||| * vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv * iattem=0 50 continue * call shell(pri,pri(lanelg),pri(lanelc),pri(lanode),pri(lacoeg), x pri(laitri),pri(laited),nelr,nrtri,nno) * call dwmatg(pri(laited),3,nrtri,0,'i','itedge') * stop * if(iattem.gt.0) call dimpr(pri,nd2,1,8,lad2ph,-1) nd2=nrtri call dimpr(pri,nd2,1,8,lad2ph,1) iattem=iattem+1 * ********************************************************* * solve ********************************************************* * call rhslhs( i pri,pri(lanelg),pri(lacoeg),pri(lanoco),pri(laitri), i pri(laited),pri(lasegd),pri(laisgd),pri(larhd), i pri(lacond),pri(laconv),pri(laabso),refchr, o pri(lau),pri(lad2ph),ul2,uene,umax,cstab) * * tolt=tol*ul2 for l2-control * tolt=tol*uene for l2-control of gradient * tolt=tol*umax for maxnorm control * if(refchr.eq.'ENE') ulevel=uene if(refchr.eq.'MAX') ulevel=umax if(refchr.eq.'RMS') ulevel=ul2 ulevel=dmax1(ulevel,1.d-10) tolt=tol*ulevel tolreft=tolref*ulevel * ********************************************************* * predict/check and modify mesh if necessary ********************************************************* * nrerr=nelr call dimpr(pri,nrerr,1,8,laerr,1) * ********************************************************* * meshsize control ********************************************************* * call mark(pri(lanode),pri(lacoeg),pri(lanelc),pri(lau), $ pri(laerr),pri(lad2ph),pri(laleve),pri(lanelg), $ pri(lairef),refchr,pri(lacond),pri(laconv), $ pri(laabso),pri(larhd),pri(laisgd),pri(lasegd), $ cstab,pri(lanoco),pri(laitri),pri(laited),espac) * ********************************************************* * Write the error to output ********************************************************* * * call test(pri(laitri),pri(lacoeg),pri(laerr),espac,nrtri) call wrerr(pri(laitri),pri(lacoeg),pri(laited),pri(lanoco), & pri(lau),espac,ulevel,isla) * * old: if(espac.lt.tolt.or.lvl.eq.maxl) then * * write(*,*) 'espac=',espac,' tolreft=',tolreft if(espac.lt.tolreft) then irefin=0 else call errcal(nrtri,pri(lanelg),pri(laerr),pri(lairef), & pri(laleve),lvl,nelr,pri(lacoeg),pri(lanode),ulevel, & mact,irefin,espac,tol,tolt,tolreft,refchr) endif * call dwmatg(pri(lairef),1,nrtri,0,'i','iref1') * call dimpr(pri,nrerr,1,8,laerr,-1) * ********************************************************* * Plot ********************************************************* * imesh=imesh+1 call dimpr(pri,nno,2,8,laarr,1) call meshpl(pri,pri(lacoeg),pri(lau), $ pri(laitri),pri(lacond),pri(laited), $ pri(lanelg),pri(laarr),imesh) call dimpr(pri,nno,2,8,laarr,-1) * ********************************************************* * if tolerance not fulfilled, change mesh ********************************************************* * if(irefin.ge.1) then * ********************************************************* * Adaptive refinement if marker irefin.ge.1 ********************************************************* * lvlbef=lvl * call arefin( $ pri(lanode),pri(lacoeg),pri(lanelc), $ pri(laleve),pri(lanelg),pri(lairef), $ pri(lanoco),pri(lavnor),pri(lau)) * lvlaft=lvl if(lvlbef.eq.lvlaft) goto 60 * * nrtri=nrtri * goto 50 endif * * ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ * ||||||||||||||||||||||||||||||||||||||||||||||||||| * 60 continue call wrsol(77,78,pri(lacoeg),pri(lau),pri(laitri), & pri(laited),imesh) * return end * * * subroutine shell(pri,nelgrp,nelcon,nodes,coeg,itri, x itedge,nelr,nrtri,nno) implicit double precision (a-h,o-z) dimension coeg(2,*),nelgrp(5,*),itri(4,*),nodes(3,*) dimension nelcon(3,*),iv(4),itedge(3,*),pri(*) * ********************************************************* * triangle shell itri ********************************************************* * nrtri=0 * do 100 iel=1,nelr if(nelgrp(1,iel).ne.0) goto 100 * do 20 i=1,3 iv(i)=nodes(i,iel) 20 continue if(nelgrp(2,iel).eq.0) then nrtri=nrtri+1 itri(1,nrtri)=iv(1) itri(2,nrtri)=iv(2) itri(3,nrtri)=iv(3) itri(4,nrtri)=iel else do 40 i=1,3 nei=nelcon(i,iel) if(nei.le.0) goto 40 ison=nelgrp(1,nei) if(ison.gt.0) then iside=i do 30 j=1,3 if(nelcon(j,nei).eq.iel) inod=nodes(j,ison) 30 continue endif 40 continue nrtri=nrtri+1 itri(4,nrtri)=iel itri(1,nrtri)=iv(iside) k=(11*iside-3*iside*iside-4)/2 itri(2,nrtri)=iv(k) itri(3,nrtri)=inod nrtri=nrtri+1 itri(4,nrtri)=iel itri(1,nrtri)=iv(iside) itri(2,nrtri)=inod k=(3*iside*iside-13*iside+16)/2 itri(3,nrtri)=iv(k) * endif 100 continue * ********************************************************* * create edge connections ********************************************************* * call dimpr(pri,4,nno*4,4,laa,1) call nodlis(nno,nrtri,pri(laa),itri,nlist,4) call edges(nodes,nelcon,nrtri,pri(laa),itri,itedge) call dimpr(pri,4,nno*4,4,laa,-1) * return end * * * subroutine nodlis(nno,nele,nodel,itri,nlist,inr) * implicit double precision (a-h,o-z) dimension nodel(4,*),itri(inr,*) * imax=nno+1 nlist=nno do 50 iel=1,nele do 40 i=1,3 inod=itri(i,iel) 10 continue ik=nodel(4,inod) if(ik.lt.0) inod=-nodel(4,inod) if(ik.lt.0) goto 10 if(ik.eq.3) goto 20 nodel(ik+1,inod)=iel nodel(4,inod)=ik+1 goto 40 20 continue nodel(4,inod)=-imax nodel(1,imax)=iel nodel(4,imax)=1 nlist=imax imax=imax+1 40 continue 50 continue return end * * * subroutine edges(nodes,nelcon,nele,nodel,itri,itedge) * implicit double precision (a-h,o-z) dimension nodes(3,*),nelcon(3,*) dimension nodel(4,*),itri(4,*),itedge(3,*) * ********************************************************* * create edge connections ********************************************************* * do 50 iel=1,nele iv1=itri(1,iel) iv2=itri(2,iel) iv3=itri(3,iel) ielfa=itri(4,iel) do 40 i=1,3 * ********************************************************* * first loop for iv2 ********************************************************* * inod1=iv2 10 continue ik1=nodel(4,inod1) inr1=ik1 if(ik1.lt.0) inr1=3 do 30 ites1=1,inr1 iel1=nodel(ites1,inod1) if(iel1.eq.iel) goto 30 * ********************************************************* * second loop for iv3 ********************************************************* * inod2=iv3 20 continue ik2=nodel(4,inod2) inr2=ik2 if(ik2.lt.0) inr2=3 do 25 ites2=1,inr2 iel2=nodel(ites2,inod2) if(iel1.eq.iel2) goto 35 25 continue if(ik2.lt.0) inod2=-ik2 if(ik2.lt.0) goto 20 30 continue if(ik1.lt.0) inod1=-ik1 if(ik1.lt.0) goto 10 do 31 k=1,3 ivt1=(11*k-3*k*k-4)/2 ivt2=(3*k*k-13*k+16)/2 ivt3=(3*k*k-13*k+16)/2 if(nodes(ivt1,ielfa).eq.iv2.and. & nodes(ivt2,ielfa).eq.iv3) iedg=k if(nodes(ivt1,ielfa).eq.iv3.and. & nodes(ivt2,ielfa).eq.iv2) iedg=k 31 continue itedge(i,iel)=nelcon(iedg,ielfa) goto 36 35 continue itedge(i,iel)=iel1 36 ii=iv1 iv1=iv2 iv2=iv3 iv3=ii 40 continue 50 continue return end * * * subroutine dimpr(pri, & n,m, & nbyte, & la, & mode) * ********************************************************* * get space for array in pri * * n number of rows * m number of columns * * nbyte number of bytes in each element * * la startposition in pri * * mode = 1 allocate space in pri * 0 print out used space in pri * -1 release arrays in pri * -2 release all arrays in pri * 2 initiate database ********************************************************* * implicit double precision (a-h,o-z) dimension pri(*) common /dbasco/ ifree,ilast * ********************************************************* * calculate number of bytes for this array ********************************************************* * iimode = mode + 3 go to ( 40,10,50,10,70),iimode 10 continue mb=n*m*nbyte if ( mod(mb,8) .ne. 0 ) mb=mb+8-mod(mb,8) mb=mb/8 if ( mode .ne. 1 ) go to 20 la=ifree ifree=ifree+mb if ( ifree-1 .gt. ilast ) go to 60 call zerop(pri(la),mb*8) go to 30 20 continue ifree=ifree-mb 30 continue * go to 80 * 40 continue ifree=1 go to 80 * 50 continue iff=ifree-1 write (6,60000) iff write (6,60010) ilast go to 80 * ********************************************************* * error ********************************************************* * 60 continue write(6,*) 'memory too small' stop 70 continue ifree=1 ilast=n * 80 continue * return 60000 format('0 last used adress in pri is ',i9) 60010 format(' allocated length is ',i9) end * * * subroutine bari(x,y,coeg,iv,c,da) implicit double precision (a-h,o-z) dimension coeg(2,*),c(3),iv(3) * ********************************************************* * bari computes the barycentric * coordinates of the point (x,y) relative to * the triangle with vertices specified in iv. * the coordinates are returned in c ********************************************************* * iv1=iv(1) iv2=iv(2) iv3=iv(3) x2=coeg(1,iv2)-coeg(1,iv1) y2=coeg(2,iv2)-coeg(2,iv1) x3=coeg(1,iv3)-coeg(1,iv1) y3=coeg(2,iv3)-coeg(2,iv1) xr=x-coeg(1,iv1) yr=y-coeg(2,iv1) det=x2*y3-x3*y2 c(2)=(xr*y3-x3*yr)/det c(3)=(x2*yr-xr*y2)/det c(1)=1.0d0-c(2)-c(3) da=dabs(det)/2.d0 return end * * * subroutine wrerr(itri,coeg,itedge,nocon,u, & errtot,ulevel,isla) * ********************************************************* * writes on file 6 : true errors (true solution * given in function exact(x,y,t) ********************************************************* * implicit double precision (a-h,o-z) * dimension u(*),itri(4,*),coeg(2,*),itedge(3,*),nocon(*) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * if(isla.eq.1) goto 10 isla=1 * write(6,100) write(6,*) 10 continue * relat=100.d0*errtot/ulevel write(6,200) lvl,nvact,nrtri,errtot,int(relat) return * * 100 format(' ','level',2x,'nodes',6x,'elements',4x, & 'approximative err.',4x,'relative app. err. (%)') 200 format(' ',i2,3x,i6,7x,i6,8x,e12.5,12x,i3) end * * * subroutine errcal(nrtri,nelgrp,err,iref,level, $ lvl,nelr,coeg,nodes,ulevel,mact,irefin,espac, $ tol,tolt,tolreft,refchr) * ********************************************************* * Calculate error ********************************************************* * implicit double precision (a-h,o-z) character*3 refchr dimension nelgrp(5,*),err(*),iref(*),level(*) dimension coeg(2,*),nodes(3,*),xc(3),yc(3) * irefin=0 * eps=1.d-4 * itesr=0 rele=nrtri/2.d0 * 10 continue errtest=0.d0 do 600 iel=1,nelr if(nelgrp(1,iel).ne.0) goto 600 * do i=1,3 xc(i)=coeg(1,nodes(i,iel)) yc(i)=coeg(2,nodes(i,iel)) enddo hold=dsqrt(2.d0*tarea(xc,yc)) errtest=errtest+err(iel)*hold*hold*tarea(xc,yc) iref(iel)=0 * ********************************************************* * err1 is jumps plus load ********************************************************* * res=dmax1(err(iel),eps) * * Maximum norm * if(refchr.eq.'MAX') then * * h^2 * Log(1/h) * Res = tolt * hnew=dsqrt((tolt/dabs(dlog(hold)))/res) endif * * L_2-norm * if(refchr.eq.'RMS') then * * h^2 * Res * area * nele = tolt * hnew=dsqrt(dsqrt(tolt/res/rele)) endif * * Energy norm * if(refchr.eq.'ENE') then * * h^2 * Res^2 * area * nele = tolt^2 * hnew=dsqrt(tolt/res/dsqrt(rele)) endif * nn=0 if(hnew.lt.hold) then nn=1 endif if(hnew.lt.hold.and.level(iel).eq.lvl) then itesr=itesr+1 endif if(hnew.gt.2.d0*hold) then nn=0 endif * iref(iel)=nn 600 continue if(itesr.gt.0) irefin=1 * * CHeck percentage * ip=INT(100*itesr/rele) * ********************************************************* * Check whether the mesh has been targeted for * refinement ********************************************************* * * write(6,*) 'irefin=',irefin * write(6,*) 'tolt=',tolt,' tolreft',tolreft * write(6,*) '**********************************' * if(irefin.eq.0.and.tolt.gt.tolreft) then if(ip.lt.5.and.tolt.gt.tolreft) then tolt=tolt/2.d0 tol=tol/2.d0 goto 10 endif tolt=tolt/2.d0 tol=tol/2.d0 * return end * * * subroutine divid2(iel1,nodes,coeg,nelcon,level,nelgrp, $ iref,nocon,u,iel2) implicit double precision (a-h,o-z) * ********************************************************* * refine element iel1 if possible ********************************************************* * dimension nodes(3,*),coeg(2,*),nelcon(3,*), $ level(*),u(*),nelgrp(5,*),iref(*), $ nocon(*),ieln(4),vn(2) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * ********************************************************* * iel1 to be refined, check connections ********************************************************* * do 10 i=1,3 neighb=nelcon(i,iel1) if (neighb.le.0) go to 10 * ********************************************************* * has neighbor, check level ********************************************************* * if (level(neighb).lt.level(iel1)) go to 80 10 continue * ********************************************************* * ok to divide ********************************************************* * iel2=iel1 do 20 i=1,4 ielfr1=ielfr ieln(i)=ielfr1 nelgrp(i,iel1)=ieln(i) do 25 j=2,4 25 nelgrp(j,ieln(i))=0 nelgrp(5,ieln(i))=iel1 level(ieln(i))=level(iel1)+1 lvl=max0(lvl,level(ieln(i))) call elfree(ielfr,nelgrp,nelr) * 20 continue * do 70 i=1,3 * ********************************************************* * fill in nodes ********************************************************* * nodes(i,ieln(i+1))=nodes(i,iel1) nelcon(i,ieln(i+1))=ieln(1) nelcon(i,ieln(1))=ieln(i+1) * neighb=nelcon(i,iel1) if (neighb.le.0) go to 30 if (nelgrp(1,neighb).eq.0) go to 30 go to 60 30 continue * ********************************************************* * neighbouring to the boundary or larger element, * new node 'inofr1' ********************************************************* * nvact=nvact+1 inofr1=inofr inofr=nocon(inofr1) if(inofr.le.nno) goto 40 nocon(inofr)=inofr+1 40 continue if(inofr1.gt.nno) nno=nno+1 if (neighb.gt.0) nelgrp(2,neighb)=nelgrp(2,neighb)+1 nodes(i,ieln(1))=inofr1 j=i*i-5*i+7 k=4*i-i*i nodes(j,ieln(k))=inofr1 j=4*i-i*i-1 k=i*i-5*i+8 nodes(j,ieln(k))=inofr1 * j=(11*i-3*i*i-2)/2 k=(3*i*i-13*i+18)/2 * * iela=ieln(j) * ielb=ieln(k) * nvface=nvface+1 * nelcon(i,iela)=neighb * nvface=nvface+1 * nelcon(i,ielb)=neighb * nelcon(i,ieln(j))=nelcon(i,iel1) nelcon(i,ieln(k))=nelcon(i,iel1) * * ********************************************************* * fill in new coordinate ********************************************************* * j=(i*i-5*i+8)/2 k=(3*i-i*i+4)/2 l=6-k-j i1=nodes(j,iel1) i2=nodes(k,iel1) * nocon(inofr1)=min0(nocon(i1),nocon(i2)) * if (neighb.gt.0) nocon(inofr1)=-neighb*10 * u(inofr1)=(u(i1)+u(i2))/2.d0 * * if(neighb.lt.0) then * call findspl( * $ 0.5d0,coeg(1,i1),coeg(1,i2), * $ vnorm(1,neighb),vnorm(3,neighb), * $ coeg(1,inofr1),vn) * else * coeg(1,inofr1)=(coeg(1,i1)+coeg(1,i2))/2.d0 coeg(2,inofr1)=(coeg(2,i1)+coeg(2,i2))/2.d0 * endif go to 70 60 continue * ********************************************************* * neighbor has the same level ********************************************************* * do 61 ij=1,3 if(nelcon(ij,neighb).eq.iel1) then inod=nodes(ij,nelgrp(1,neighb)) goto 62 endif 61 continue 62 continue * nocon(inod)=0 nodes(i,ieln(1))=inod j=i*i-5*i+7 k=4*i-i*i nodes(j,ieln(k))=inod j=4*i-i*i-1 k=i*i-5*i+8 nodes(j,ieln(k))=inod * j=(11*i-3*i*i-2)/2 k=(3*i*i-13*i+18)/2 * do 64 ij=1,3 if(nelcon(ij,neighb).eq.iel1) then l=(11*ij-3*ij*ij-2)/2 m=(3*ij*ij-13*ij+18)/2 in1=nelgrp(l,neighb) in2=nelgrp(m,neighb) nelcon(i,ieln(k))=in1 nelcon(i,ieln(j))=in2 do 63 jj=1,3 if(nelcon(jj,in1).eq.iel1) nelcon(jj,in1)=ieln(k) if(nelcon(jj,in2).eq.iel1) nelcon(jj,in2)=ieln(j) 63 continue goto 65 endif 64 continue 65 continue * 70 continue go to 90 80 continue * ********************************************************* * not possible to divide ********************************************************* * iel2=neighb return * 90 continue ir=iref(iel1)-1 iref(iel1)=0 do 100 i=1,4 100 iref(nelgrp(i,iel1))=ir return end * * * subroutine mark(nodes,coeg,nelcon,u,err,d2ph, & level,nelgrp,iref,refchr,cond,conv,abso,rhs, & isegd,segdat,cstab,nocon,itri,itedge,espac) * ********************************************************* * mark triangles to be refined ********************************************************* * implicit double precision (a-h,o-z) character*3 refchr dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*), $ err(*),nelgrp(5,*),iref(*),nocon(*),u(*),d2ph(*), $ itri(4,*),itedge(3,*),cond(*),conv(*),abso(*), $ rhs(*),isegd(*),segdat(2,*) * common /comnav/ $ tol,digits, $ mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, $ inofr,nrgre,nvact,nrtri,nelact common /poisson/ nrhs,ncond,nconv,nabso,numseg * error=0.d0 * do 10 iel=1,nrtri ielfa=macrotr(iel,itri,nelgrp) * call jumps(iel,ielfa,emax,u,itri,coeg,nocon, $ cond,ncond,conv,nconv,abso,nabso,isegd,segdat, $ refchr,d2ph,cstab,rhs,nrhs,itedge,res,error) * ielf=itri(4,iel) err(ielf)=dmax1(res,err(ielf)) 10 continue * call dwmatg(err,1,nrtri,0,'rd','err') * if(refchr.eq.'ENE') then espac=dsqrt(error) else espac=error endif * return end * * * subroutine refine(iel,nodes,coeg,nelcon, & level,nelgrp,iref,nocon,u) implicit double precision (a-h,o-z) dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*), & nelgrp(5,*),iref(*),nocon(*),u(*) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * ********************************************************* * iel to be refined, check if possible ********************************************************* * if(level(iel).ge.maxl) then iref(iel)=0 return endif * 10 continue iel1=iel call divide( $ iel1,nodes,coeg,nelcon,level, $ nelgrp,iref,nocon,u,iel2) * if (iel1.eq.iel2) go to 30 20 continue iel1=iel2 call divide( $ iel1,nodes,coeg,nelcon,level, $ nelgrp,iref,nocon,u,iel2) * if (iel2.ne.iel1) go to 20 go to 10 30 continue return end * * * subroutine elfree(ielfr,nelgrp,nelr) implicit double precision (a-h,o-z) * ********************************************************* * find new free element number and change nelgrp ********************************************************* * dimension nelgrp(5,*) * ielfr1=ielfr ielfr=iabs(nelgrp(1,ielfr)) if(ielfr.gt.nelr) nelgrp(1,ielfr)=-(ielfr+1) if(ielfr1.gt.nelr) nelr=nelr+1 nelgrp(1,ielfr1)=0 return end * * * subroutine arefin(nodes,coeg,nelcon, & level,nelgrp,iref,nocon,vnorm,u) implicit double precision (a-h,o-z) * ********************************************************* * refine according to iref(iel) ********************************************************* * dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*), $ nelgrp(5,*),iref(*),nocon(*),vnorm(4,*),u(*) * common /comnav/ $ tol,digits, $ mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, $ inofr,nrgre,nvact,nrtri,nelact * istart=1 imax=nelr do 40 iel=istart,imax * if(iref(iel).gt.0) call refine(iel,nodes,coeg,nelcon, & level,nelgrp,iref,nocon,u) 40 continue * do 30 iel=1,nelr if (iref(iel).lt.0) call reduce(iel,nodes,coeg,nelcon, & level,nelgrp,iref,nocon,u) 30 continue * ********************************************************* * fix green nodes ********************************************************* * call grdiv(nodes,coeg,nelcon,u,nocon, & level,nelgrp,iref) return end * * * subroutine reduce(iel1,nodes,coeg,nelcon, & level,nelgrp,iref,nocon,u) implicit double precision (a-h,o-z) * ********************************************************* * delete group containing element iel1 if possible ********************************************************* * dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*), & nelgrp(5,*),iref(*),nocon(*),u(*),ieln(4) common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * ********************************************************* * iel1 to be unrefined, check * connections of the group ********************************************************* * ifa=nelgrp(5,iel1) if(ifa.eq.0) goto 90 do 20 ie=1,4 iels=nelgrp(ie,ifa) if(iref(iels).gt.0) goto 90 ieln(ie)=iels do 10 i=1,3 neighb=nelcon(i,iels) if (neighb.le.0) go to 10 * ********************************************************* * has neighbor, check level ********************************************************* * if (nelgrp(1,neighb).ne.0) go to 90 10 continue 20 continue * ********************************************************* * ok to reduce ********************************************************* * imid=nelgrp(1,ifa) do 12 ij=1,4 12 nelgrp(ij,ifa)=0 do 13 ij=1,3 neifa=nelcon(ij,ifa) if(neifa.le.0) goto 13 if(nelgrp(1,neifa).gt.0) nelgrp(2,ifa)=nelgrp(2,ifa)+1 if(nelgrp(1,neifa).eq.0) nelgrp(2,neifa)=nelgrp(2,neifa)-1 13 continue do 30 i=1,4 nelgrp(1,ieln(i))=-ielfr ielfr=ieln(i) do 30 j=2,4 nelgrp(j,ieln(i))=0 30 continue * ********************************************************* * test for neighbors to father ********************************************************* * do 70 ie=1,3 inod=nodes(ie,imid) neighb=nelcon(ie,ifa) if (neighb.le.0) go to 40 if (nelgrp(1,neighb).eq.0) go to 40 go to 50 40 continue * ********************************************************* * neighbouring to the boundary or larger element, * delete node ********************************************************* * nvact=nvact-1 nocon(inod)=inofr u(inod)=0.d0 inofr=inod go to 70 50 continue * ********************************************************* * neighbor has the same level, node becomes green ********************************************************* * nocon(inod)=-ifa*10 * ********************************************************* * change connections of neighbours ********************************************************* * do 55 ij=1,3 if(nelcon(ij,neighb).eq.ifa) then l=(11*ij-3*ij*ij-2)/2 m=(3*ij*ij-13*ij+18)/2 nelcon(ij,nelgrp(l,neighb))=ifa nelcon(ij,nelgrp(m,neighb))=ifa goto 56 endif 55 continue 56 continue 60 continue 70 continue * ********************************************************* * blank out nodes array ********************************************************* * do 80 i=1,4 do 80 j=1,3 80 nodes(j,ieln(i))=0 * go to 100 90 continue * ********************************************************* * not possible to unrefine ********************************************************* * iref(iel1)=0 return * 100 continue do 110 i=1,4 110 iref(ieln(i))=0 return end * * * subroutine calch(iel,nodes,coeg,h) implicit double precision (a-h,o-z) * ********************************************************* * compute local meshsize h ********************************************************* * dimension nodes(3,*),coeg(2,*) h=0.d0 do 10 i=1,2 i1=i+1 x1=coeg(1,nodes(i,iel)) y1=coeg(2,nodes(i,iel)) x2=coeg(1,nodes(i1,iel)) y2=coeg(2,nodes(i1,iel)) dl=dsqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)) h=dmax1(dl,h) 10 continue h=h/1.414213d0 return end * * * function tarea(xc,yc) implicit double precision (a-h,o-z) * ********************************************************* * compute area of triangle with vertex * coordinates in (xc,yc) ********************************************************* * dimension xc(*),yc(*),dl(3) * do 10 i=1,3 j=(11*i-3*i*i-4)/2 dx=xc(j)-xc(i) dy=yc(j)-yc(i) dl(i)=dsqrt(dx*dx+dy*dy) 10 continue * s=(dl(1)+dl(2)+dl(3))/2.d0 tarea=dsqrt(s*(s-dl(1))*(s-dl(2))*(s-dl(3))) return end * * * subroutine uiter(u,uite,nv,eps) implicit double precision (a-h,o-z) * dimension u(*),uite(*) * alfa=0.d0 beta=0.d0 do 10 i=1,nv alfa=alfa+dabs(u(i)) beta=beta+dabs(u(i)-uite(i)) 10 continue * eps=beta/alfa * return end * * * subroutine knots(iel,iv,itri) implicit double precision (a-h,o-z) * dimension itri(4,*),iv(*) * iv(1)=itri(1,iel) iv(2)=itri(2,iel) iv(3)=itri(3,iel) * return end * * * subroutine zerop(vec,ntal) dimension vec(*) n=ntal/4 do 10 i=1,n 10 vec(i)=0. return end * * * subroutine mvcp(rout,rin,ntal) dimension rout(*),rin(*) n=ntal/4 do 10 i=1,n 10 rout(i)=rin(i) return end * * * subroutine grad(ux,uy,coeg,u,iv) implicit double precision (a-h,o-z) * dimension coeg(2,*),u(*),iv(*) * ********************************************************* * this routine computes the gradient of the * scalar grid function u in the triangle with * nodes given in iv ********************************************************* * iv1=iv(1) iv2=iv(2) iv3=iv(3) x2=coeg(1,iv2)-coeg(1,iv1) x3=coeg(1,iv3)-coeg(1,iv1) y2=coeg(2,iv2)-coeg(2,iv1) y3=coeg(2,iv3)-coeg(2,iv1) u2=u(iv2)-u(iv1) u3=u(iv3)-u(iv1) det=x2*y3-x3*y2 ux=(u2*y3-u3*y2)/det uy=(x2*u3-x3*u2)/det return end * * * subroutine jumps(iel,ielfa,emax,u,itri,coeg,nocon, & cond,ncond,conv,nconv,abso,nabso,isegd,segdat, & refchr,d2ph,cstab,rhs,nrhs,itedge,res,error) * ********************************************************* * compute difference quotients in element iel ********************************************************* * implicit double precision (a-h,o-z) character*3 refchr dimension itri(4,*),coeg(2,*),itedge(3,*), & u(*),nocon(*),cond(ncond,*),rhs(nrhs,*),d2ph(*), & conv(nconv,*),abso(nabso,*),isegd(*),segdat(2,*) dimension iv(3),xc(3),yc(3),fi(3),fix(3),fiy(3) dimension condu(2,2),veloc(2) character*80 rhscha,cndcha(3),cnvcha(2),abscha,ex * common /chardata/ idrhs,idcnd,idcnv,idabs, $ rhscha,cndcha,cnvcha,abscha * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * ********************************************************* * convective residual ********************************************************* * xm=0.d0 ym=0.d0 um=0.d0 reserr=0.d0 resjum=0.d0 ux=0.d0 uy=0.d0 * do 10 ijk=1,3 iv(ijk)=itri(ijk,iel) xc(ijk)=coeg(1,iv(ijk)) yc(ijk)=coeg(2,iv(ijk)) xm=xm+coeg(1,iv(ijk))/3.d0 ym=ym+coeg(2,iv(ijk))/3.d0 um=um+u(iv(ijk))/3.d0 10 continue * call base(xm,ym,xc,yc,fi,fix,fiy,area) hloc=dsqrt(2.d0*area) * * ********************************************************* * absorption ********************************************************* * if(idabs.eq.0) call evalexpr(abscha, xm, ym, prod, ierr) if(idabs.eq.1) call shepardv( i xm,ym,nabso,abso(1,1),abso(1,2),abso(1,3),1, o prod) if(idabs.eq.2) prod=abso(1,ielfa) * ********************************************************* * convection ********************************************************* * if(idcnv.eq.0) then ex=cnvcha(1) call evalexpr(ex, xm, ym, veloc(1), ierr) ex=cnvcha(2) call evalexpr(ex, xm, ym, veloc(2), ierr) endif if(idcnv.eq.1) call shepardv( i xm,ym,nconv,conv(1,1),conv(1,2),conv(1,3),2, o veloc) if(idcnv.eq.2) then veloc(1)=conv(1,ielfa) veloc(2)=conv(2,ielfa) endif * ********************************************************* * load ********************************************************* * if(idrhs.eq.0) call evalexpr(rhscha, xm, ym, fxy, ierr) if(idrhs.eq.1) call shepardv( i xm,ym,nrhs,rhs(1,1),rhs(1,2),rhs(1,3),1, o fxy) if(idrhs.eq.2) fxy=rhs(1,ielfa) * velx=veloc(1) vely=veloc(2) do 20 i=1,3 reserr=reserr+(velx*fix(i)+vely*fiy(i))*u(iv(i)) reserr=reserr+prod*fi(i)*u(iv(i)) ux=ux+fix(i)*u(iv(i)) uy=uy+fiy(i)*u(iv(i)) 20 continue reserr=dabs(reserr-fxy) * ********************************************************* * calculate influence from laplacian ********************************************************* * iv1=iv(1) iv2=iv(2) iv3=iv(3) * uxb=ux uyb=uy do 90 in=1,3 call gaussx(3,in,xc,yc,xm,ym) itnei=itedge(in,iel) if(itnei.lt.0) isname= -itnei if(isegd(isname).gt.0) then gload=segdat(1,isname) gamma=segdat(2,isname) endif * ********************************************************* * dirichlet boundary -> skip this side ********************************************************* * if(itnei.le.0.and.isegd(isname).lt.0) goto 89 * ********************************************************* * conduction ********************************************************* * if(idcnd.eq.0) then ex=cndcha(1) call evalexpr(ex, xm, ym, condu(1,1), ierr) ex=cndcha(2) call evalexpr(ex, xm, ym, condu(1,2), ierr) ex=cndcha(3) call evalexpr(ex, xm, ym, condu(2,2), ierr) condu(2,1)=condu(1,2) endif if(idcnd.eq.1) then call shepardv( i xm,ym,ncond,cond(1,1),cond(1,2),cond(1,3),3, o condu(1,1)) condu(2,2)=condu(1,2) condu(1,2)=condu(2,1) endif if(idcnd.eq.2) then condu(1,1)=cond(1,ielfa) condu(2,1)=cond(2,ielfa) condu(1,2)=condu(2,1) condu(2,2)=cond(3,ielfa) endif * ux1 = condu(1,1)*uxb+condu(1,2)*uyb uy = condu(2,1)*uxb+condu(2,2)*uyb ux = ux1 * r1=-(coeg(2,iv3)-coeg(2,iv2)) r2=coeg(1,iv3)-coeg(1,iv2) h=sqrt(r1*r1+r2*r2) r1=r1/h r2=r2/h * ********************************************************* * solution at neighbor ********************************************************* * un=r1*ux+r2*uy unn=0.d0 if(itnei.gt.0) then call grad(uxn,uyn,coeg,u,itri(1,itnei)) ux1 = condu(1,1)*uxn+condu(1,2)*uyn uyn = condu(2,1)*uxn+condu(2,2)*uyn uxn = ux1 unn=r1*uxn+r2*uyn else unn=-(gload-gamma*um) endif * ********************************************************* * compute jump between iel and itnei * resjum is ********************************************************* * ed1=0.5d0*(un-unn)/hloc resjum=resjum+ed1*ed1*h * 89 ii=iv1 iv1=iv2 iv2=iv3 iv3=ii 90 continue resjum=dsqrt(resjum*hloc/area) * ********************************************************* * interpolation constants and total residual ********************************************************* * c1=2.d0/7.d0 c2=2.d0/5.d0 res=c1*reserr+c2*resjum * if(refchr.eq.'MAX') then res=res/10.d0 temp=cstab*hloc*hloc*res error=dmax1(error,temp) error=error res=res*cstab endif * if(refchr.eq.'RMS') then argu=hloc*hloc*res*d2ph(iel) error=error+argu*area res=res*d2ph(iel) endif * if(refchr.eq.'ENE') then res=res/10.d0 resol=res det=condu(1,1)*condu(2,2)-condu(1,2)*condu(2,1) condu(1,2) = -condu(1,2)/det condu(2,1) = -condu(2,1)/det temp = condu(1,1)/det condu(1,1) = condu(2,2)/det condu(2,2) = temp res1 = (condu(1,1)+condu(1,2))*res/2.d0 res2 = (condu(2,1)+condu(2,2))*res/2.d0 res = dsqrt(res1*res1+res2*res2) error=error+(hloc*hloc*resol*res)*area endif return end * * * subroutine grdiv(nodes,coeg,nelcon,u,nocon, & level,nelgrp,iref) implicit double precision (a-h,o-z) * * clear up too many green knots * dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*), & nocon(*),nelgrp(5,*),iref(*),u(*) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * 100 continue do 10 iel=1,nelr if(nelgrp(1,iel).ne.0) goto 10 if(nelgrp(2,iel).eq.0) goto 10 * nrgr=nelgrp(2,iel) * ********************************************************* * if more than one green node, divide ********************************************************* * if(nrgr.lt.2) goto 10 call refine(iel,nodes,coeg,nelcon, & level,nelgrp,iref,nocon,u) goto 100 10 continue return end subroutine test(itri,coeg,err,espac,nrtri) implicit double precision (a-h,o-z) dimension coeg(2,*),err(*),itri(4,*),xc(3),yc(3) errtot=0.d0 do 10 iel=1,nrtri do 5 i=1,3 xc(i)=coeg(1,itri(i,iel)) yc(i)=coeg(2,itri(i,iel)) 5 continue area=tarea(xc,yc) h=dsqrt(2.d0*area) res=err(iel) * errtot=errtot+h*h*res*res*area temp=h*h*res*dabs(dlog(h)) errtot=dmax1(errtot,temp*temp) 10 continue * write(6,*) 'espac=',espac * write(6,*) 'errtot=',dsqrt(errtot) return end subroutine fooread(file,vnorm,itedge,nface) implicit none double precision vnorm(4,*) integer file,itedge(4,*),nele,nno,nface,i,j * read(file,*) nele,nno,nface read(file,*) ((itedge(i,j),i=1,3),j=1,nele) read(file,*) ((vnorm(i,j),i=1,4),j=1,nface) close(unit=file) return end subroutine divide(iel1,nodes,coeg,nelcon,level,nelgrp, & iref,nocon,u,iel2) implicit double precision (a-h,o-z) * ********************************************************* * refine element iel1 if possible ********************************************************* * dimension nodes(3,*),coeg(2,*),nelcon(3,*),level(*),u(*), & nelgrp(5,*),iref(*),nocon(*),ieln(4) * common /comnav/ & tol,digits, & mact,nelr,nno,ielfr,neltot,maxt,maxv,maxl,lvl, & inofr,nrgre,nvact,nrtri,nelact * ********************************************************* * iel1 to be refined, check connections ********************************************************* * do 10 i=1,3 neighb=nelcon(i,iel1) if (neighb.le.0) go to 10 * ********************************************************* * has neighbor, check level ********************************************************* * if (level(neighb).lt.level(iel1)) go to 80 10 continue * ********************************************************* * ok to divide ********************************************************* * iel2=iel1 do 20 i=1,4 ielfr1=ielfr ieln(i)=ielfr1 nelgrp(i,iel1)=ieln(i) do 25 j=2,4 25 nelgrp(j,ieln(i))=0 nelgrp(5,ieln(i))=iel1 level(ieln(i))=level(iel1)+1 lvl=max0(lvl,level(ieln(i))) call elfree(ielfr,nelgrp,nelr) * 20 continue * do 70 i=1,3 * ********************************************************* * fill in nodes ********************************************************* * nodes(i,ieln(i+1))=nodes(i,iel1) nelcon(i,ieln(i+1))=ieln(1) nelcon(i,ieln(1))=ieln(i+1) * neighb=nelcon(i,iel1) if (neighb.le.0) go to 30 if (nelgrp(1,neighb).eq.0) go to 30 go to 60 30 continue * ********************************************************* * neighbouring to the boundary or larger element, * new node 'inofr1' ********************************************************* * nvact=nvact+1 inofr1=inofr inofr=nocon(inofr1) if(inofr.le.nno) goto 40 nocon(inofr)=inofr+1 40 continue if(inofr1.gt.nno) nno=nno+1 if (neighb.gt.0) nelgrp(2,neighb)=nelgrp(2,neighb)+1 nodes(i,ieln(1))=inofr1 j=i*i-5*i+7 k=4*i-i*i nodes(j,ieln(k))=inofr1 j=4*i-i*i-1 k=i*i-5*i+8 nodes(j,ieln(k))=inofr1 * j=(11*i-3*i*i-2)/2 k=(3*i*i-13*i+18)/2 nelcon(i,ieln(j))=nelcon(i,iel1) nelcon(i,ieln(k))=nelcon(i,iel1) * * ********************************************************* * fill in new coordinate ********************************************************* * j=(i*i-5*i+8)/2 k=(3*i-i*i+4)/2 l=6-k-j nocon(inofr1)=min0(nocon(nodes(j,iel1)),nocon(nodes(k,iel1))) * if (neighb.gt.0) nocon(inofr1)=-neighb*10 * u(inofr1)=0.5d0*(u(nodes(j,iel1))+u(nodes(k,iel1))) * do 50 ij=1,2 coeg(ij,inofr1)=0.5d0*(coeg(ij,nodes(j,iel1))+ & coeg(ij,nodes(k,iel1))) 50 continue go to 70 60 continue * ********************************************************* * neighbor has the same level ********************************************************* * do 61 ij=1,3 if(nelcon(ij,neighb).eq.iel1) then inod=nodes(ij,nelgrp(1,neighb)) goto 62 endif 61 continue 62 continue * nocon(inod)=0 nodes(i,ieln(1))=inod j=i*i-5*i+7 k=4*i-i*i nodes(j,ieln(k))=inod j=4*i-i*i-1 k=i*i-5*i+8 nodes(j,ieln(k))=inod * j=(11*i-3*i*i-2)/2 k=(3*i*i-13*i+18)/2 * do 64 ij=1,3 if(nelcon(ij,neighb).eq.iel1) then l=(11*ij-3*ij*ij-2)/2 m=(3*ij*ij-13*ij+18)/2 in1=nelgrp(l,neighb) in2=nelgrp(m,neighb) nelcon(i,ieln(k))=in1 nelcon(i,ieln(j))=in2 do 63 jj=1,3 if(nelcon(jj,in1).eq.iel1) nelcon(jj,in1)=ieln(k) if(nelcon(jj,in2).eq.iel1) nelcon(jj,in2)=ieln(j) 63 continue goto 65 endif 64 continue 65 continue * 70 continue go to 90 80 continue * ********************************************************* * not possible to divide ********************************************************* * iel2=neighb return * 90 continue ir=iref(iel1)-1 iref(iel1)=0 do 100 i=1,4 100 iref(nelgrp(i,iel1))=ir return end