# to unbundle, sh this file (in an empty directory) mkdir src echo src/mgghat.f 1>&2 sed >src/mgghat.f <<'//GO.SYSIN DD src/mgghat.f' 's/^-//' - subroutine mgghat - include 'commons' - real nwctrg,mgf0,tguess - external clock -c -c perform initializations -c - mgf0 = mgfreq - call init -c -c do refine/solve loop until relative error estimate is small enough -c - 1 if (rerest .lt. tol .or. ierr.ne.0) go to 2 -c -c timing of this step -c - timetl = clock(0) -c -c refinement/solution/error estimation step -c - call refine - call solve - call esterr -c -c finish timing of this step -c - timetl = clock(0) - timetl - if (outlev.ge.2) write(ioutpt,100) timetl,clock(0)-timett -c -c estimate the amount of time needed for another refine/solve step, -c and if there isn't enough time left, modify mgfreq or quit -c - tguess = timetl * mgfreq - if (tguess + clock(0)-timett .gt. mxtime) then - if (timetl .gt. 10.*r1mach(4)) then - nwctrg = (mxtime - (clock(0)-timett))/timetl - else - nwctrg = mgfreq - endif - if (nwctrg.le.1.) then - ierr=5 - else - if (nwctrg.lt.mgfreq) mgfreq=nwctrg - endif - endif - go to 1 -c -c end of main loop -c - 2 continue - mgfreq = mgf0 -c -c ending output -c - call outend -c - 100 format(/' time for this refinement/solution step',f10.2/ - 1 ' total time so far ',f10.2) -c -c continue to view the plots until the menu quit button is pressed -c - if (grafic) then - if (menuon) then - write(6,1100) - grpaws = .false. - 11 if (.not.menuon) go to 12 - call mn2mg - go to 11 - 12 continue - else - write(6,1200) - read * - endif -c -c clear all the graphics in case the program is called again -c - menuon = .false. - grafic = .false. - do 21 i=1,100 - if (pltsel(i)) then - pltsel(i)=.false. - if (gunit(i).ne.-1) then - call gpclos(gunit(i)) - gunit(i)=-1 - endif - endif - 21 continue - endif - 1100 format(/'The Graphics Selections Menu is still active, and can'/ - . 'be used to add, delete or rotate gnuplot displays.'// - . 'To terminate the program, press Quit on the graphics menu.') - 1200 format(/' The Graphics will remain on screen until you', - . ' press return.') -c - return - end -c -c block data to initialize parameters to default values. The -c user can override the defaults with assignment statements in -c the main program. -c - block data - include 'commons' -c -c Stopping critera related to space allocation. Default is the -c parameters used to dimension arrays -c - data mxvert,mxtri,mxlev,mxnode / ndvert,ndtri,ndlev,ndnode / -c -c Maximum allowed cpu time in seconds, as a stopping criteria. -c Default is 12 hours -c - data mxtime / 43200. / -c -c Error tolerance, as a stopping criteria. Returns if the relative -c energy norm error estimate drops below tol. Default is 0., i.e., -c use some other stopping criteria. -c - data tol / 0. / -c -c I/O unit to which to write printed output. Default is 6. -c - data ioutpt / 6 / -c -c I/O unit for gnuplot file output. Default is 4. -c - data gpfile / 4 / -c -c Amount of printed output. Usable values are: -c <= 0 No output, except error messages -c 1 Header plus summary at end of execution -c 2 Information after each phase of execution -c 3 Low level of debug information -c 4 Medium level of debug information -c >= 5 High level of debug information -c Values of 3, 4 and 5 are probably not useful to normal users. -c Default is 2. -c - data outlev / 2 / -c -c Order of the polynomial approximation (degree+1). -c 2 is linear, 3 is quadratic and 4 is cubic. Higher orders are -c possible, but require code modification for more accurate -c quadrature rules. Default is 4. -c - data iorder / 4 / -c -c parameters to pin down a nonunique solution -c nuniq = .true. if the user supplies the value of the solution at some -c vertex -c nuniqx, nuniqy = the x,y coordinates where the solution is given -c nuniqv = the value of the solution at that vertex -c If the user says it's not unique, but doesn't give the point and/or -c value, the default is to use the value 0. at the point (0.,0.) -c - data nuniq / .false. / - data nuniqx, nuniqy, nuniqv / 0.,0.,0. / -c -c Uniform (.true.) or adaptive (.false.) refinement. Default is .false. -c - data unifrm / .false. / -c -c Flags for which gnuplot files to write. Default is none. -c gptri = 0 no triangulation file -c 1 write triangulation file -c gpsol = 0 no solution file -c = n write file with solution on n X n grid -c gpconv= 0 no convergence file -c 1 write convergence file -c - data gptri, gpsol, gpconv / 0, 0, 0 / -c -c do/do_not prompt the user during initialization for selection -c of graphics or startup of graphics menu. Default is to prompt. -c - data gquiet / .false. / -c -c Use a widget based menu for graphics selection. Default is false. -c - data menuon / .false. / -c -c pltsel is a logical array of plot selections. The i'th entry is true -c if graphical displays are to be made for the corresponding entry -c in the following table. Default is all false. -c 1 Computed Solution; Surface 31 Error; Surface -c 2 Computed Solution; Contour 32 Error; Contour -c 3 Computed Solution; Facets 33 Error; Facets -c 4 Computed Solution; Surf & Tri 34 Error; Surf & Tri -c 6 Computed Solution; Facet & Tri 36 Error; Facet & Tri -c 11 True Solution; Surface 41 Triangulation -c 12 True Solution; Contour 51 Nodes vs. energy err -c 13 True Solution; Facets 52 Nodes vs. max error -c 14 True Solution; Surf & Tri 53 Nodes vs. err est -c 16 True Solution; Facet & Tri 54 Nodes vs. err & est -c 21 Comp & True Solutions; Surface 61 Time vs. energy err -c 22 Comp & True Solutions; Contour 62 Time vs. max error -c 23 Comp & True Solutions; Facets 63 Time vs. err est -c 24 Comp & True Solutions; Surf & Tri 64 Time vs. err & est -c 26 Comp & True Solutions; Facet & Tri -c - data pltsel / 100*.false. / -c -c The remaining parameters determine details of the numerical -c method used. Modifying these would only be of interest to -c a numerical analyst experimenting with adaptive multilevel -c methods. -c -c How much to multiply the number of nodes by in each refinement phase, -c i.e., amount of refinement between multigrid phases. -c - data mgfreq / 2. / -c -c Number of relaxation iterations before and after coarse grid correction -c - data nu1, nu2 / 1, 2 / -c -c number of multigrid cycles in each solution phase -c - data ncyc / 1 / -c - end -c -c -------- init -c - subroutine init - include 'commons' - integer i,holdol - logical holdgr - external clock -c -c header output -c - if (outlev.ge.1) call outhed -c -c set up graphics -c - if (.not. gquiet) call setupg -c -c initial triangulation -c - call inittr -c -c check validity of input -c - call valinp -c -c initializations -c - if (outlev.ge.2) write(ioutpt,100) - 100 format(/' begin initialization') -c -c start clock for total time -c - timett = clock(0) - timert = 0. - timest = 0. - timeet = 0. -c -c set the list of vertices of each level to be empty -c and the number of vertices of each level to 0 -c - do 10 i=2,mxlev - frstvt(i) = -1 - lvert (i) = 0 - lbvert(i) = 0 - 10 continue -c -c initially 1 level and no error -c - nvert0=nvert - nlev=1 - ierr=0 - gplev=1 -c -c initialize basis values and quadrature formula -c - call initbs -c -c set node renumbering tables -c - call initrn -c -c initialize triangulation data structures -c - call initds -c -c initialize nodes -c - call initnd -c -c set up equations for coarsest grid -c - call initeq -c -c initialize for exact coarse grid solution -c - call initxs -c -c solve coarse grid problem -c - call exsolv - if (outlev.ge.4) then - write(ioutpt,300) - call prvec(u,1,1) - endif - 300 format(/' initial solution'/) -c -c compute initial errors -c - call unorm(un,trun) - unrm = un - trunrm = trun - call errors(emaxv,emaxn,emaxq,eenrg) - gerr = eenrg - if (trunrm.ge.1e-10) then - rgerr = gerr/trunrm - else - rgerr = gerr - endif - emax = emaxq - if (emaxn.gt.emax) emax=emaxn -c -c save information for gnuplot convergence file -c - gpmxer(gplev)=emax - gpener(gplev)=rgerr - gpnode(gplev)=nnode -c -c compute first error indicators -c - holdol = outlev - if (outlev.le.2) outlev = 0 - holdgr = grafic - grafic = .false. - call esterr - outlev = holdol - grafic = holdgr -c -c initialize graphics -c - call initgr -c -c compute time used for initialization -c - timei = clock(0) - timett - if (outlev.ge.2) write(ioutpt,200) timei - 200 format(/' initializations complete'// - 1 ' time for initialization ',f10.2) -c - return - end -c -c -------- initbs -c - subroutine initbs - include 'commons' - integer i,j,sumord,sub,inc,nod,deg,nodebc(3),qp,zeta,rnod - real rdeg,nf,factor -c -c initializations associated with basis functions -c - if (outlev.ge.3) write(ioutpt,400) - 400 format(' initializations associated with bases') -c -c number of new nodes associated with a vertex (interior and boundary) -c and number of nodes in a triangle -c - nnodev=(iorder-1)**2 - nnodvb=(iorder*(iorder-1))/2 - nnodtr=(iorder*(iorder+1))/2 -c -c quadrature rule -c - if (iorder.eq.2) then - nqpt=1 - quadpt(1,1)=1./3. - quadpt(2,1)=1./3. - quadpt(3,1)=1./3. - quadw(1)=1. - nqptb=1 - qptb(1)=.5 - qwtb(1)=1. - elseif (iorder.eq.3) then - nqpt=3 - quadpt(1,1)=2./3. - quadpt(2,1)=1./6. - quadpt(3,1)=quadpt(2,1) - quadw(1)=1./3. - quadpt(1,2)=quadpt(2,1) - quadpt(2,2)=quadpt(1,1) - quadpt(3,2)=quadpt(2,1) - quadw(2)=quadw(1) - quadpt(1,3)=quadpt(2,1) - quadpt(2,3)=quadpt(2,1) - quadpt(3,3)=quadpt(1,1) - quadw(3)=quadw(1) - nqptb=2 - qptb(1)=.5+1./(2.*sqrt(3.)) - qptb(2)=1.-qptb(1) - qwtb(1)=.5 - qwtb(2)=.5 - elseif (iorder.eq.4) then - nqpt=6 - quadpt(1,1)=.8168475729 - quadpt(2,1)=.0915762135 - quadpt(3,1)=quadpt(2,1) - quadw(1)=.1099517436 - quadpt(1,2)=quadpt(2,1) - quadpt(2,2)=quadpt(1,1) - quadpt(3,2)=quadpt(2,1) - quadw(2)=quadw(1) - quadpt(1,3)=quadpt(2,1) - quadpt(2,3)=quadpt(2,1) - quadpt(3,3)=quadpt(1,1) - quadw(3)=quadw(1) - quadpt(1,4)=.1081030181 - quadpt(2,4)=.4459484909 - quadpt(3,4)=quadpt(2,4) - quadw(4)=.2233815896 - quadpt(1,5)=quadpt(2,4) - quadpt(2,5)=quadpt(1,4) - quadpt(3,5)=quadpt(2,4) - quadw(5)=quadw(4) - quadpt(1,6)=quadpt(2,4) - quadpt(2,6)=quadpt(2,4) - quadpt(3,6)=quadpt(1,4) - quadw(6)=quadw(4) - nqptb=3 - qptb(1)=(1.+sqrt(.6))/2. - qptb(2)=.5 - qptb(3)=1.-qptb(1) - qwtb(1)=5./18. - qwtb(2)=4./9. - qwtb(3)=qwtb(1) - elseif (nqpt.le.0) then - write(ioutpt,100) iorder - stop - 100 format(/' ********FATAL ERROR********'// - 1 ' need a quadrature rule which is compatable with order ',i2) - endif -c -c number of basis changes associated with a node -c - sumord=(iorder*(iorder+1))/2 - do 10 i=1,nnodev - nbasch(i)=sumord - 10 continue - sub=1 - inc=iorder-1 - 1 nbasch(sub)=iorder - sub=sub+inc - inc=inc-2 - if (inc.gt.0) go to 1 - inc=iorder-2 - if (iorder.eq.2) go to 3 - 2 nbasch(sub)=iorder - sub=sub+inc - inc=inc-2 - if (inc.gt.0) go to 2 - 3 continue -c -c set index of old nodes for new nodes for basis change -c - do 20 i=1,nnodvb - do 20 j=1,nbasch(i) - ibasch(j,i)=j - 20 continue - if (iorder.eq.2) go to 51 - do 50 i=nnodvb+1,nnodev - do 30 j=1,iorder - ibasch(j,i)=j - 30 continue - do 40 j=iorder+1,sumord - ibasch(j,i)=sumord-iorder+j - 40 continue - 50 continue - 51 continue -c -c set barycentric coordinates of red nodes in old triangles (scaled -c by 2*degree) -c note this is only used here and olndvt is used as a temporary for this -c - deg=iorder-1 - rdeg=float(deg) - nod=1 - do 60 i=1,iorder/2 - do 60 j=2*deg+1-2*i,2*i-1,-2 - olndvt(nod,1)=j - olndvt(nod,2)=2*i-1 - olndvt(nod,3)=2*deg-olndvt(nod,1)-olndvt(nod,2) - nod=nod+1 - 60 continue - if (iorder.eq.2) go to 71 - do 70 i=1,deg/2 - do 70 j=2*deg+1-2*i,2*i+1,-2 - olndvt(nod,1)=2*i-1 - olndvt(nod,2)=j - olndvt(nod,3)=2*deg-olndvt(nod,1)-olndvt(nod,2) - nod=nod+1 - 70 continue - 71 continue -c -c evaluate basis function -c -c need value and derivatives at the quadrature points for quadratures -c and the values at red nodes for basis changes -c - nodebc(1)=iorder - nodebc(2)=-1 - nodebc(3)=0 -c -c for each node of a triangle -c - do 240 nod=1,nnodtr -c -c find barycentric coordinates of node (scaled by degree) -c - nodebc(1)=nodebc(1)-1 - nodebc(2)=nodebc(2)+1 - if (nodebc(1).lt.0) then - nodebc(3)=nodebc(3)+1 - nodebc(2)=0 - nodebc(1)=deg-nodebc(3) - endif -c -c initialize basis values at 1. and derivatives at 0. -c -c normalization factor - nf=1. -c quadrature points - do 120 qp=1,nqpt - qpbas(nod,qp)=1. - do 110 zeta=1,3 - qpdbdz(zeta,nod,qp)=0. - 110 continue - 120 continue -c red nodes - do 130 rnod=1,nnodvb - cbasch(nod,rnod)=1. - 130 continue -c -c evaluate basis -c - do 180 i=1,3 - do 170 j=0,deg-1 - if (nodebc(i).gt.j) then -c adjust normalization factor - nf=nf*float(nodebc(i)-j)/rdeg -c values and derivatives at quadrature points - do 150 qp=1,nqpt -c new factor in basis value - factor=quadpt(i,qp)-float(j)/rdeg -c derivatives with respect to barycentric coordinates - do 140 zeta=1,3 - qpdbdz(zeta,nod,qp)=qpdbdz(zeta,nod,qp)*factor - 140 continue - qpdbdz(i,nod,qp)=qpdbdz(i,nod,qp)+qpbas(nod,qp) -c basis value - qpbas(nod,qp)=qpbas(nod,qp)*factor - 150 continue -c value at red nodes - do 160 rnod=1,nnodvb - cbasch(nod,rnod)=cbasch(nod,rnod)* - 1 float(olndvt(rnod,i)-2*j)/(2.*rdeg) - 160 continue - endif - 170 continue - 180 continue -c -c boundary quadrature points -c - if (nod.le.iorder) then - do 181 qp=1,nqptb - qpbasb(nod,qp)=1. - qpbbdz(nod,qp)=0. -c basis value - do 171 j=0,deg - if (j.ne.nod-1) then - qpbasb(nod,qp)=qpbasb(nod,qp)* - & (1.-qptb(qp)-float(j)/rdeg)/(float(nod-1-j)/rdeg) - endif - 171 continue -c derivative - do 161 j=0,deg - if (j.ne.nod-1) then - qpbbdz(nod,qp)=qpbbdz(nod,qp)+qpbasb(nod,qp) - & /(1.-qptb(qp)-float(j)/rdeg) - endif - 161 continue - 181 continue - endif -c -c normalize values -c - do 210 qp=1,nqpt - qpbas(nod,qp)=qpbas(nod,qp)/nf - do 190 zeta=1,3 - qpdbdz(zeta,nod,qp)=qpdbdz(zeta,nod,qp)/nf - 190 continue - 210 continue - do 220 rnod=1,nnodvb - cbasch(nod,rnod)=cbasch(nod,rnod)/nf - 220 continue -c -c copy values at red nodes for matching triangle -c - if (iorder.eq.2) go to 231 - sub=1 - do 230 rnod=nnodvb+1,nnodev - 229 sub=sub+1 - if (nbasch(sub).eq.iorder) go to 229 - cbasch(nod,rnod)=cbasch(nod,sub) - 230 continue - 231 continue - 240 continue -c -c compress 0's out of basis change values -c @future do this more carefully. I have checked that for -c iorder=2,3 and 4 all the 0's are at the end, so this works. -c don't know if it works for higher order -c may also want to try to not put the 0's there in the first place -c - do 310 i=1,nnodev - j=nbasch(i) - 311 if (cbasch(j,i).ne.0.) go to 312 - j=j-1 - go to 311 - 312 continue - nbasch(i)=j - 310 continue -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,500) - 500 format(/' basis initialization complete') - write(ioutpt,501) nnodev,nnodvb,nnodtr - 501 format(' number of new nodes associated with a vertex'/ - 1 ' interior ',i3,' boundary ',i3/ - 2 ' number of nodes in a triange ',i3) - write(ioutpt,502) - 502 format(' number of basis changes for each red node') - write(ioutpt,503) (nbasch(i),i=1,nnodev) - 503 format(1x,19i4) - write(ioutpt,504) - 504 format(' barycentric coordinates of red nodes') - write(ioutpt,505) ((olndvt(i,j),j=1,3),i=1,nnodvb) - 505 format(1x,3i5) - write(ioutpt,506) - 506 format(' index and constant for basis changes') - do 510 i=1,nnodev - write(ioutpt,511) i - 511 format(' red node ',i3) - write(ioutpt,512) (ibasch(j,i),cbasch(j,i),j=1,nbasch(i)) - 512 format(1x,4(i4,1pe15.8)) - 510 continue - write(ioutpt,513) - 513 format(' quadrature points and weights') - write(ioutpt,514) ((quadpt(j,i),j=1,3),quadw(i),i=1,nqpt) - 514 format(1x,4(1pe15.8)) - write(ioutpt,533) - 533 format(' boundary quadrature points and weights') - write(ioutpt,534) (qptb(i),qwtb(i),i=1,nqptb) - 534 format(1x,2(1pe15.8)) - write(ioutpt,515) - 515 format(' values and derivatives of bases at quadrature points') - do 520 i=1,nnodtr - write(ioutpt,521) i - 521 format(' basis at node ',i3) - do 520 j=1,nqpt - write(ioutpt,522) qpbas(i,j),qpdbdz(1,i,j),qpdbdz(2,i,j), - 1 qpdbdz(3,i,j) - 522 format (1x,4(1pe15.8)) - 520 continue - write(ioutpt,535) - 535 format(' values of bases at boundary quadrature points') - do 540 i=1,iorder - write(ioutpt,536) i,(qpbasb(i,j),j=1,nqptb) - 536 format(' node ',i3,2x,4(1pe15.8)) - 540 continue - write(ioutpt,635) - 635 format(' values of derivatives at boundary quadrature points') - do 640 i=1,iorder - write(ioutpt,536) i,(qpbbdz(i,j),j=1,nqptb) - 640 continue - endif - if(outlev.ge.3) write(ioutpt,650) - 650 format(' basis initializations complete') - return - end -c -c -------- initds -c - subroutine initds - include 'commons' - integer i,j,k,vert -c -c initialize triangulation data structures -c - if (outlev.ge.3) write(ioutpt,100) - 100 format(' initialize triangulation data structures') -c -c neighbors of each triangle -c - do 110 i=1,ntri - k=1 - 111 if (k.gt.ntri) go to 112 - if(vertex(2,k).eq.vertex(2,i) .and. - 1 vertex(3,k).eq.vertex(3,i) .and. i.ne.k) go to 112 - k=k+1 - go to 111 - 112 continue - if (k.le.ntri) then - neigh(1,i)=k - else - if (neigh(1,i).ge.0) then - write(ioutpt,200) 1,i - stop - endif - endif - k=1 - 113 if (k.gt.ntri) go to 114 - if(vertex(1,k).eq.vertex(1,i) .and. - 1 vertex(3,k).eq.vertex(3,i) .and. i.ne.k) go to 114 - k=k+1 - go to 113 - 114 continue - if (k.le.ntri) then - neigh(2,i)=k - else - if (neigh(2,i).ge.0) then - write(ioutpt,200) 2,i - stop - endif - endif - k=1 - 115 if (k.gt.ntri) go to 116 - if(vertex(1,k).eq.vertex(1,i) .and. - 1 vertex(2,k).eq.vertex(2,i) .and. i.ne.k) go to 116 - k=k+1 - go to 115 - 116 continue - if (k.le.ntri) then - neigh(3,i)=k - else - if (neigh(3,i).ge.0) then - write(ioutpt,200) 3,i - stop - endif - endif - 110 continue - 200 format(' ********FATAL ERROR********'/ - . ' side opposite vertex ',i1,' of triangle ',i6/ - . ' is on the boundary but does not have neigh assigned', - . ' the boundary piece.') -c -c triangles around each vertex -c - do 10 i=1,nvert - do 10 j=1,8 - tringl(j,i)=0 - 10 continue - do 20 i=1,ntri - do 20 j=1,3 - vert=vertex(j,i) - k=1 - 1 if (tringl(k,vert).eq.0) go to 2 - k=k+1 - if (k.gt.8) then - write(ioutpt,101) vert - stop - 101 format(/' ********FATAL ERROR********'// - 1 ' too many triangles around vertex ',i6) - endif - go to 1 - 2 continue - tringl(k,vert)=i - 20 continue -c -c determine which vertices are boundary -c boundary vertices occur in triangles with a nonpositive neighbor -c - do 30 i=1,nvert - vrtlev(i)=1 - 30 continue - do 40 i=1,ntri - if (neigh(1,i).le.0) then - vrtlev(vertex(2,i))=-1 - vrtlev(vertex(3,i))=-1 - endif - if (neigh(2,i).le.0) then - vrtlev(vertex(1,i))=-1 - vrtlev(vertex(3,i))=-1 - endif - if (neigh(3,i).le.0) then - vrtlev(vertex(1,i))=-1 - vrtlev(vertex(2,i))=-1 - endif - 40 continue -c -c set up linked lists of vertices -c and count number of vertices -c - lvert(1)=0 - lbvert(1)=0 - frstvt(1)=-1 - do 50 i=1,nvert - nextvt(i)=frstvt(1) - frstvt(1)=i - if (vrtlev(i).gt.0) then - lvert(1)=lvert(1)+1 - else - lbvert(1)=lbvert(1)+1 - endif - 50 continue -c -c debug output -c - if (outlev.ge.4) then - call plttri(1,1) - endif - if (outlev.ge.3) write(ioutpt,201) - 201 format(' triangulation data structures initialization complete') -c - return - end -c -c -------- initeq -c - subroutine initeq - include 'commons' - integer i,j,t,c,r -c -c initialize equations, i.e., matrix and right side -c and set up initial error indicator problems -c - if (outlev.ge.3) write(ioutpt,100) - 100 format(' initialize equations') -c -c set maximum number of nonzeroes in lower and upper -c triangular parts of a row. mxidup is sum of lower and upper -c - mxidlo = 2*iorder*(iorder-1) - mxidup = 3*mxidlo -c -c set coef and rs by doing integrals over each -c triangle and assembling -c -c initialize to 0 -c - do 120 i=1,nnode - rs(i)=0. - do 110 j=1,mxidlo - coef(j,i)=0. - idcoef(j,i)=0 - 110 continue - coef(mxidlo+1,i)=0. - do 115 j=mxidlo+1,mxidup - idcoef(j,i)=0 - 115 continue - 120 continue -c - do 140 t=1,ntri -c quadratures for this triangle - call quad(t,1) -c make sure column pointers are all there - do 125 i=1,nadd - if (row(i).eq.col(i)) go to 125 - c=1 - 121 if (c.gt.mxidlo .or. idcoef(c,row(i)).eq.col(i) - 1 .or. idcoef(c,row(i)).eq.0) go to 122 - c=c+1 - go to 121 - 122 continue - if (c.gt.mxidlo) then - write(ioutpt,500) - stop - 500 format(' ********FATAL ERROR********'// - 1 ' ran out of room in lower idcoef in initeq') - else - idcoef(c,row(i))=col(i) - endif - 125 continue -c add values to matrix - call addmat -c add values to right side - do 130 i=1,naddrs - rs(rowrs(i))=rs(rowrs(i))+addrs(i) - 130 continue - 140 continue -c -c set upper triangular column pointers -c - do 150 i=1,nnode - do 150 j=1,mxidlo - if (idcoef(j,i).ne.0) then - r=idcoef(j,i) - c=mxidlo+1 - 141 if (c.gt.mxidup .or. idcoef(c,r).eq.i .or. - 1 idcoef(c,r).eq.0) go to 142 - c=c+1 - go to 141 - 142 continue - if (c.gt.mxidup) then - write(ioutpt,501) - stop - 501 format(' ********FATAL ERROR********'// - 1 ' ran out of room in upper idcoef in initeq') - else - idcoef(c,r)=i - endif - endif - 150 continue -c -c set up initial error indicator problems -c -c set lists of triangles to be empty -c - do 209 i=1,4 - eihead(i)=-1 - eitail(i)=-1 - 209 continue -c -c set up equations -c - do 210 i=1,ntri - if (neigh(3,i).le.0 .or. neigh(3,i).gt.i) then - call setei(i) - else - call setei2(neigh(3,i),i) - endif - 210 continue -c -c debug output -c - if (outlev.ge.4) call outmat - if (outlev.ge.3) write(ioutpt,511) - 511 format(' matrix initialization complete') -c - return - end -c -c -------- initgr -c - subroutine initgr - include 'commons' -c -c initializations for graphics -c -c set all gnuplot units to -1 -c - do 15 i=1,100 - gunit(i)=-1 - 15 continue -c -c set grafic to true iff one of pltsel is true, or menuon is true -c - grafic = .false. - if (menuon) grafic = .true. - do 20 i=1,100 - if (pltsel(i)) grafic = .true. - 20 continue -c -c set graphics pause (grpaws) true, to allow time for -c reshaping windows -c - grpaws = .true. -c -c if runtime graphics are on, and the number of isolines for -c surface plots was not specified by the user, set it to 20 -c -c set number of isolines for surface plots to default 20X20 -c - gpsolx = 20 - gpsoly = 20 -c -c create data files for starting plots, check for new message from -c the menu, and draw the plots -c - if (grafic) then - call filtri - call filsol - call filcon - if (menuon) call mn2mg - call grftri - call grfsol - call grfcon - endif -c - return - end -c -c ---- setupg -c - subroutine setupg - include 'commons' -c -c set up the run time graphics -c - logical there - character*1 inchar - character*200 inline - character*10 digits - integer ipoint, select, inval, inlen - data digits / '0123456789' / -c - write(ioutpt,100) - 100 format(' ') - write(ioutpt,210) - 210 format(/' NOTE -- you must have gnuplot installed and in'/ - . ' your execution path for graphics to work'/) - write(ioutpt,110) - 110 format('Do you want any runtime graphics (y/n)?') - read(5,'(a1)') inchar - if (inchar .eq. 'y') then - grafic = .true. - write(ioutpt,220) - 220 format(/' NOTE -- you must have tcl/tk installed and wish'/ - . ' in your execution path for the menu to ', - . 'work'/) - write(ioutpt,120) - 120 format('Do you want menu-widget driven selection (y/n)?') - read(5,'(a1)') inchar - if (inchar .eq. 'y') then - inquire(file=tmpdir//'men2mgg',exist=there) - if (there) call system('rm '//tmpdir//'men2mgg') - call system('wish < grmenu &') - menuon = .true. - write(ioutpt,130) - 130 format(/ - . ' A process has been spawned to run the widget based' - . ,' menu.'/' From the menu you can add and delete displays,'/ - . ' rotate the view of 3D plots, and change the number of'/ - . ' isolines used for surface plots. You can continue to '/ - . ' use the menu as the program executes, and after the '/ - . ' completion of the solution process. The initial '/ - . ' selections made now will not appear until the program '/ - . ' continues, and the displays may be slow to update.'// - . ' >>> Make initial selections from the menu and press ', - . 'return.'/) - read * - else - write(ioutpt,140) -c -c initialize for read -c - read(5,'(a200)') inline - inlen = len(inline) - ipoint=0 - select = 0 -c -c loop while in an integer -c - 1 continue - ipoint = ipoint + 1 - if (ipoint.gt.inlen) go to 4 - read(inline(ipoint:ipoint),'(a1)',end=4) inchar - inval = index(digits,inchar)-1 - if (inval .lt. 0 .or. inval .gt. 9) go to 2 - select = 10*select+inval - go to 1 -c -c end of integer; set pltsel -c - 2 continue - if (select .gt. 0 .and. select .lt. 100) - . pltsel(select) = .true. -c -c loop while not in integer -c - 3 continue - ipoint = ipoint + 1 - if (ipoint.gt.inlen) go to 4 - read(inline(ipoint:ipoint),'(a1)',end=4) inchar - select = index(digits,inchar)-1 - if (select .ge. 0 .and. select .le. 9) go to 1 - go to 3 -c -c end of input string -c - 4 continue - endif - endif - 140 format(/' Enter a list of integers to select any number of '/ - . ' displays from the following list:'// - . ' 1 Computed Solution; Surface 31 Error; Surface'/ - . ' 2 Computed Solution; Contour 32 Error; Contour'/ - . ' 3 Computed Solution; Facets 33 Error; Facets'/ - . ' 4 Computed Solution; Surf & Tri 34 Error; Surf & Tri'/ - . ' 6 Computed Solution; Facet & Tri 36 Error; Facet & Tri'/ - . '11 True Solution; Surface 41 Triangulation'/ - . '12 True Solution; Contour 51 Nodes vs. energy err'/ - . '13 True Solution; Facets 52 Nodes vs. max error'/ - . '14 True Solution; Surf & Tri 53 Nodes vs. err est'/ - . '16 True Solution; Facet & Tri 54 Nodes vs. err & est'/ - . '21 Comp & True Solutions; Surface 61 Time vs. energy err'/ - . '22 Comp & True Solutions; Contour 62 Time vs. max error'/ - . '23 Comp & True Solutions; Facets 63 Time vs. err est'/ - . '24 Comp & True Solutions; Surf & Tri 64 Time vs. err & est'/ - . '26 Comp & True Solutions; Facet & Tri'/) -c - return - end -c -c -------- initnd -c - subroutine initnd - include 'commons' - integer t,nay,i,j,sub,inc,k,deg,lnod,v1,v2,nod,bctype - real x,y - real bcrhs,cu -c -c define nodes for the initial triangulation -c - if (outlev.ge.3) write(ioutpt,100) - 100 format(' initialize nodes') -c - nnode = nvert - do 60 t=1,ntri -c -c nodes which are vertices -c - node(1,t) = vertex(1,t) - node(iorder,t) = vertex(2,t) - node(nnodtr,t) = vertex(3,t) - if (iorder.gt.2) then -c -c nodes along base -c - nay = neigh(3,t) - if (nay.gt.t .or. nay.le.0) then - do 10 i=2,iorder-1 - nnode = nnode+1 - node(i,t) = nnode - if (nay.gt.0) node(i,nay) = nnode - 10 continue - endif -c -c nodes along older side -c - nay = neigh(2,t) - if (nay.gt.t .or. nay.le.0) then - sub = iorder+1 - inc = iorder-1 - do 20 i=1,iorder-2 - nnode = nnode + 1 - node(sub,t) = nnode - if (nay.gt.0) node(sub,nay) = nnode - sub = sub + inc - inc = inc - 1 - 20 continue - endif -c -c nodes along newer side -c - nay = neigh(1,t) - if (nay.gt.t .or. nay.le.0) then - sub = 2*iorder-1 - inc = iorder-2 - do 30 i=1,iorder-2 - nnode = nnode + 1 - node(sub,t) = nnode - if (nay.gt.0) node(sub,nay) = nnode - sub = sub + inc - inc = inc - 1 - 30 continue - endif -c -c nodes interior to triangle -c - if (iorder.gt.3) then - sub = iorder + 2 - do 50 i=1,iorder-3 - do 40 j=1,iorder-2-i - nnode = nnode + 1 - node(sub,t) = nnode - sub = sub + 1 - 40 continue - sub = sub + 2 - 50 continue - endif - endif - 60 continue -c -c set number of initial nodes -c - nnode0 = nnode -c -c set boundary node flags and boundary condition -c - do 110 i=1,nnode - bndnod(i) = .false. - u(i) = 0. - 110 continue -c -c if the solution is not unique (due to pure Neumann or periodic -c conditions), then set the single equation at a vertex whose -c coordinates were given by the user to tie down the solution -c - if (nuniq) then -c -c find the vertex -c - eps = 4.*r1mach(4) - do 115 i=1,nvert - if (abs(xvert(i)-nuniqx).lt.eps .and. - . abs(yvert(i)-nuniqy).lt.eps) then - ii = i - go to 116 - endif - 115 continue -c -c didn't find a vertex with the right coordinates -c - write(ioutpt,200) nuniqx,nuniqy - 200 format(/' ********FATAL ERROR********'// - . 'Non-unique solution must be specified at a vertex.'/ - . 'There is no vertex with coordinates',2(1pe15.8)) - if (outlev.le.2) then - write(ioutpt,201) - else - write(ioutpt,202) - write(ioutpt,203) (xvert(i),yvert(i),i=1,nvert) - endif - stop - 201 format(/'To see a list of vertex coordinates, rerun with ', - . 'outlev=3') - 202 format(/'The vertices are:'/) - 203 format(2(1pe15.8)) -c -c found the vertex, call it a Dirichlet point -c - 116 continue - bndnod(ii) = .true. - u(ii) = nuniqv - endif -c - deg = iorder-1 -c -c for each neighbor of each triangle, if the neighbor is nonpositive then -c that side of the triangle is on the boundary. mark those nodes -c - do 140 i=1,ntri -c - do 130 j=1,3 - if (neigh(j,i).le.0) then - if (j.eq.1) then - lnod= iorder - inc = deg - v1 = vertex(2,i) - v2 = vertex(3,i) - elseif (j.eq.2) then - lnod= 1 - inc = iorder - v1 = vertex(1,i) - v2 = vertex(3,i) - else - lnod= 1 - inc = 1 - v1 = vertex(1,i) - v2 = vertex(2,i) - endif -c - do 120 k=0,deg - nod = node(lnod,i) - if (.not. bndnod(nod)) then - x=((deg-k)*xvert(v1)+k*xvert(v2))/float(deg) - y=((deg-k)*yvert(v1)+k*yvert(v2))/float(deg) - call bcond(x,y,-neigh(j,i),cu,bcrhs,bctype) - if (bctype.eq.1) then - bndnod(nod)=.true. - u(nod)=bcrhs - endif - endif - lnod=lnod+inc - inc=inc-1 - if (inc.lt.1) inc=1 - 120 continue - endif - 130 continue - 140 continue -c - if (outlev.ge.4) call plttri(3,0) - if (outlev.ge.4) call outtri - if (outlev.ge.3) write(ioutpt,102) - 102 format(' node initialization complete') -c - return - end -c -c -------- initrn -c - subroutine initrn - include 'commons' - integer newnn,ordd2,i,j,sub,inc,inc2,disp,disp2,inode -c -c initialize data structures used for renumbering nodes -c -c renum contains the renumbering of the nodes ina a local numbering -c when a pair of triangles is divided. for the figures the example -c of a cubic basis is used. before division the local numbering is -c -c peak 10 -c older triangle 8 9 -c 5 6 7 -c oldest vertex 1 2 3 4 -c -c oldest vertex 1 2 3 4 -c 5 6 7 -c younger triangle 8 9 -c peak 10 -c -c when the triangles are divided, new nodes are added in the -c order (* are old nodes) -c -c peak * -c older triangle * 3 * -c * 2 * 6 * -c oldest vertex * 1 * 4 * 5 * -c * 7 * 9 * -c younger triangle * 8 * -c peak * -c -c after dividing the local numbering of the nodes is -c -c 4 4 -c 3 7 7 3 -c 2 6 9 9 6 2 -c 1 5 8 10 10 8 5 1 -c -c 1 5 8 10 10 8 5 1 -c 2 6 9 9 6 2 -c 3 7 7 3 -c 4 4 -c -c renum gives the correspondence between the new local number and -c the old local number -c -c renum(i,.) is the correspondence for the ith new triangle -c 1) the half of the older triangle with the oldest vertex (upper left) -c 2) the other half of the older triangle (upper right) -c 3) the half of the younger triangle with the oldest vertex (lower left) -c 4) the other half of the younger triangle (lower right) -c -c renum(i,j) tells how to find the jth node in the new local numbering -c for the ith triangle -c if positive, it contains the local number of the node in the -c appropriate old triangle -c if negative, it contains the negative of which new node is in -c that position -c -c with the subscripts given by the third figure above, we have for cubics -c -c 10 10 -c 8 -3 -3 8 -c 5 -2 6 6 -6 7 -c 1 -1 2 -4 -4 3 -5 4 -c -c 1 -1 2 -4 -4 3 -5 4 -c 5 -7 6 6 -9 7 -c 8 -8 -8 9 -c 10 10 -c -c - if (outlev.ge.3) write(ioutpt,100) - 100 format(' initialize node renumbering table') -c - newnn = 0 - ordd2=int(float(iorder)/2.+.6) -c -c half of the oldest triangle with the oldest node -c - inode = 1 - do 40 j=1,ordd2 - sub = j - inc = iorder - do 20 i=1,iorder-2*(j-1) - renum(1,inode) = sub - inode = inode + 1 - sub = sub + inc - inc = inc - 1 - 20 continue - if (2*j .le. iorder) then - do 30 i=1,iorder-2*j+1 - newnn = newnn - 1 - renum(1,inode) = newnn - inode = inode + 1 - 30 continue - endif - 40 continue -c -c other half of older triangle -c - disp = iorder - 1 - inc2 = iorder - 3 - inode = 1 - do 130 j=1,ordd2 - sub = iorder + 1 - j - inc = iorder - 1 - do 110 i = 1,iorder-2*(j-1) - renum(2,inode) = sub - inode = inode + 1 - sub = sub + inc - inc = inc - 1 - 110 continue - if (2*j .le. iorder) then - if (2*j .le. iorder-1) then - do 120 i=1,iorder-2*j - newnn = newnn - 1 - renum(2,inode) = newnn - inode = inode + 1 - 120 continue - endif - renum(2,inode) = - disp - inode = inode + 1 - disp = disp + inc2 - inc2 = inc2 - 2 - endif - 130 continue -c -c half of the younger triangle with the oldest node -c - disp = 1 - inc2 = iorder - 1 - inode = 1 - do 240 j=1,ordd2 - sub = j - inc = iorder - do 220 i=1,iorder-2*(j-1) - renum(3,inode) = sub - inode = inode + 1 - sub = sub + inc - inc = inc - 1 - 220 continue - if (2*j .le. iorder) then - renum(3,inode) = -disp - inode = inode + 1 - disp = disp + inc2 - inc2 = inc2 - 2 - if (2*j .le. iorder-1) then - do 230 i=1,iorder-2*j - newnn = newnn - 1 - renum(3,inode) = newnn - inode = inode + 1 - 230 continue - endif - endif - 240 continue -c -c other half of younger triangle -c - if (iorder .eq. 2*int(iorder/2)) then - disp = 1 + (iorder*iorder)/4 - else - disp = 1 + (iorder*iorder-1)/4 - endif - inc2 = iorder - 2 - disp2 = (iorder*iorder + iorder - 4)/2 - inode = 1 - do 330 j=1,ordd2 - sub = iorder + 1 - j - inc = iorder - 1 - do 310 i = 1,iorder-2*(j-1) - renum(4,inode) = sub - inode = inode + 1 - sub = sub + inc - inc = inc - 1 - 310 continue - if (2*j .lt. iorder) then - renum(4,inode) = -disp - inode = inode + 1 - disp = disp + inc2 - inc2 = inc2 - 2 - if (2*j .le. iorder-2) then - do 320 i=1,iorder-2*j-1 - newnn = newnn - 1 - renum(4,inode) = newnn - inode = inode + 1 - 320 continue - endif - renum(4,inode) = -disp2 - inode = inode + 1 - disp2 = disp2 + inc2 - endif - 330 continue - if (iorder .eq. 2*int(iorder/2)) then - renum(4,inode) = -(iorder*iorder)/4 - endif -c - if (outlev.ge.3) write(ioutpt,101) - 101 format(' renumbering table complete') -c - return - end -c -c ------- valinp -c - subroutine valinp - include 'commons' - integer t,t2,v,iv,icount -c -c check the validity of user input -c - 99 format(' ********FATAL ERROR********'/) -c -c check the order of the method -c - if (iorder.lt.2) then - write(ioutpt,99) - write(ioutpt,100) iorder - 100 format(' order of the method=',i6,' must be at least 2') - stop - endif - if (iorder.gt.4) then - write(ioutpt,200) iorder,nqpt - 200 format(' ********WARNING********'/ - 1 ' iorder=',i6,' is greater than 4.'/ - 2 ' using user provided quadrature rule with',i6,' points') - endif - if (iorder.gt.ndord) then - write(ioutpt,99) - write(ioutpt,201) iorder,ndord - 201 format(' iorder=',i6,' is .gt. ndord=',i4) - stop - endif -c -c make sure maximum allowed values do not exceed dimensions -c - if (mxvert.gt.ndvert) then - write(ioutpt,251) mxvert,ndvert - 251 format(' ********WARNING********'/ - . ' mxvert=',i10,' is greater than ndvert=',i10/ - . ' resetting mxvert equal to ndvert') - mxvert=ndvert - endif - if (mxtri.gt.ndtri) then - write(ioutpt,252) mxtri,ndtri - 252 format(' ********WARNING********'/ - . ' mxtri=',i10,' is greater than ndtri=',i10/ - . ' resetting mxtri equal to ndtri') - mxtri=ndtri - endif - if (mxnode.gt.ndnode) then - write(ioutpt,253) mxnode,ndnode - 253 format(' ********WARNING********'/ - . ' mxnode=',i10,' is greater than ndnode=',i10/ - . ' resetting mxnode equal to ndnode') - mxnode=ndnode - endif - if (mxlev.gt.ndlev) then - write(ioutpt,254) mxlev,ndlev - 254 format(' ********WARNING********'/ - . ' mxlev=',i6,' is greater than ndlev=',i6/ - . ' resetting mxlev equal to ndlev') - mxlev=ndlev - endif -c -c check that number of vertices increases during refinement -c - if (mgfreq.le.1.) then - write(ioutpt,99) - write(ioutpt,300) mgfreq - 300 format(' mgfreq=',1pe15.8,' is .le. 1.') - stop - endif -c -c check that the initial triangulation doesn't exceed dimensions -c - if (nvert.gt.ndvert .or. nvert.le.0) then - write(ioutpt,99) - write(ioutpt,351) nvert,ndvert - 351 format(' initial nvert=',i6,' is .gt. ndvert=',i6,' or .le. 0') - stop - endif - if (ntri.gt.ndtri .or. ntri.le.0) then - write(ioutpt,99) - write(ioutpt,352) ntri,ndtri - 352 format(' initial ntri=',i6,' is .gt. ndtri=',i6,' or .le. 0') - stop - endif - if ((iorder-1)*(iorder-1)*nvert.gt.ndnode) then - write(ioutpt,99) - write(ioutpt,353) (iorder-1)*(iorder-1)*nvert,ndnode - 353 format(' estimated initial nnode=',i6,' is .gt. ndnode=',i6) - stop - endif -c -c check vertex alignment and count -c - do 20 t=1,ntri - do 20 v=1,3 - iv=vertex(v,t) - if (iv.le.0 .or. iv.gt.nvert) then - write(ioutpt,99) - write(ioutpt,400) v,t,iv,nvert - 400 format(' vertex(',i6,',',i6,')=',i6, - 1 ' is .le.0 or .gt. nvert=',i6) - stop - endif - icount=0 - do 10 t2=t,ntri - if (v.eq.1) then - if (vertex(2,t2).eq.iv .or. vertex(3,t2).eq.iv) then - write(ioutpt,99) - write(ioutpt,500) iv - stop - endif - endif - if (v.eq.2) then - if (vertex(1,t2).eq.iv .or. vertex(3,t2).eq.iv) then - write(ioutpt,99) - write(ioutpt,500) iv - stop - endif - endif - if (v.eq.3) then - if (vertex(1,t2).eq.iv .or. vertex(2,t2).eq.iv) then - write(ioutpt,99) - write(ioutpt,500) iv - stop - endif - endif - 500 format(' vertex',i6,' occurs in 2 columns of vertex') - if (vertex(v,t2).eq.iv) icount=icount+1 - 10 continue - if ((v.eq.3.and.icount.gt.4).or.(v.ne.3.and.icount.gt.8)) then - write(ioutpt,99) - write(ioutpt,600) iv - 600 format(' vertex',i6,' occurs in too many triangles') - stop - endif - 20 continue -c - return - end -c -c -------- refine -c - subroutine refine - include 'commons' - integer t,limit - external clock -c -c refine triangulation -c - if (outlev.ge.2) write(ioutpt,100) - 100 format(/' begin refinement') -c -c begin timing -c - timerl = clock(0) -c -c --- uniform refinement --- -c - if (unifrm) then -c - nlev=nlev+1 - if (nlev.gt.mxlev) then - ierr=3 - nlev=nlev-1 - return - endif - limit = ntri - do 20 t=1,limit - if (iabs(vrtlev(vertex(3,t))).ne.nlev) call divtri(t) - 20 continue -c -c --- adaptive refinement --- -c - else -c -c refine until number of vertices is doubled -c - ntarg = int(mgfreq*nvert) - if (ntarg.le.nvert) ntarg=nvert+1 -c -c get triangle to divide -c - 1 call divnxt(t) -c -c divide triangle -c - call divtri(t) -c -c repeat if necessary -c - if (nvert .lt. ntarg .and. ierr .eq. 0) go to 1 -c - endif -c -c finish timing -c - timerl = clock(0) - timerl - timert = timert + timerl -c -c write data file for run time graphics, check the menu, and -c update the plots -c - if (outlev.ge.2 .and. grafic) write(ioutpt,200) - 200 format(' begin update to triangulation graphics') - if (grafic) then - call filtri - if (menuon) call mn2mg - call grftri - endif - if (outlev.ge.2 .and. grafic) write(ioutpt,300) - 300 format(' graphics complete') -c -c debug output -c - if(outlev.ge.2) call outref - if(outlev.ge.4) then - call outtri - call outmat - endif -c - return - end -c -c -------- divnxt -c - subroutine divnxt(t) - include 'commons' - integer t -c -c determine next triangle to divide -c -c if top error indicator list is empty, shift lists -c - 1 if (eihead(1) .ne. -1) go to 2 - eihead(1) = eihead(2) - eihead(2) = eihead(3) - eihead(3) = eihead(4) - eihead(4) = -1 - eitail(1) = eitail(2) - eitail(2) = eitail(3) - eitail(3) = eitail(4) - eitail(4) = -1 - if (outlev.ge.4) write(ioutpt,100) - eimax=eimax*mgfreq**(-iorder/2.) - 100 format(' shift error indicator lists for smaller eimax') - go to 1 - 2 continue -c -c return first triangle on top list -c - t = eihead(1) -c -c verify that the triangle isn't too small -c - xdiff = abs(xvert(vertex(1,t))-xvert(vertex(2,t))) - temp = abs(xvert(vertex(1,t))-xvert(vertex(3,t))) - if (temp.gt.xdiff) xdiff=temp - temp = abs(xvert(vertex(3,t))-xvert(vertex(2,t))) - if (temp.gt.xdiff) xdiff=temp - ydiff = abs(yvert(vertex(1,t))-yvert(vertex(2,t))) - temp = abs(yvert(vertex(1,t))-yvert(vertex(3,t))) - if (temp.gt.ydiff) ydiff=temp - temp = abs(yvert(vertex(3,t))-yvert(vertex(2,t))) - if (temp.gt.ydiff) ydiff=temp - xmax = abs(xvert(vertex(1,t))) - temp = abs(xvert(vertex(2,t))) - if (temp.gt.xmax) xmax=temp - temp = abs(xvert(vertex(3,t))) - if (temp.gt.xmax) xmax=temp - ymax = abs(yvert(vertex(1,t))) - temp = abs(yvert(vertex(2,t))) - if (temp.gt.ymax) ymax=temp - temp = abs(yvert(vertex(3,t))) - if (temp.gt.ymax) ymax=temp - if (xdiff/xmax.lt.10.*r1mach(4) .or. - . ydiff/ymax.lt.10.*r1mach(4)) then - call elstrm(t) - temp=errind(t) - errind(t)=0. - call elstad(t) - errind(t)=temp - go to 1 - endif -c - return - end -c -c -------- dtpair -c - subroutine dtpair(t1,t2) - include 'commons' - integer t1,t2,t(4),i,limit -c -c divide the compatably divisible pair of triangles t1,t2 -c t2<=0 if base of t1 is on boundary -c -c check for too many vertices, triangles or levels -c - if (nvert.ge.mxvert-2) then - ierr = 1 - elseif ((t2.gt.0 .and. ntri.ge.mxtri-2) - 1 .or. (t2.le.0 .and. ntri.ge.mxtri-1) ) then - ierr = 2 - elseif (iabs(vrtlev(vertex(3,t1))).ge.mxlev) then - ierr = 3 - elseif ((t2.gt.0 .and. nnode+nnodev.gt.mxnode) - 1 .or. (t2.le.0 .and. nnode+nnodvb.gt.mxnode)) then - ierr = 4 - endif - if (ierr.ne.0) return -c -c debug output -c - if (outlev.ge.4) write(ioutpt,101) t1,t2 - 101 format(' divide triangle pair ',2i6) -c -c -c begin correction of existing matrix and right side values by -c subtracting off the old inner products over the triangles divided -c - call oldip(t1,t2) -c -c adjust the data structures for vertices and triangles to divide -c - if (t2.gt.0) then - call adjds(t1,t2) - else - call adjdsb(t1) - endif -c -c finish correction fo existing matrix and right side values by -c adding in the new inner products over the new triangles -c - call newip(t1,t2) -c -c set new equation in the matrix and right side and set first -c approximate solution at the new vertex -c - call neweq(t1,t2) -c -c set up the error indicator problem for the four new triangles -c Three cases -- boundary, -c old (problem already defined in neighbor triangle) and normal. -c - if (t2.gt.0) then - t(1) = t1 - t(2) = t2 - t(3) = ntri-1 - t(4) = ntri - limit= 4 - else - t(1) = t1 - t(2) = ntri - limit= 2 - endif - do 10 i=1,limit - if (neigh(3,t(i)).le.0) then - call setei(t(i)) - elseif (neigh(3,neigh(3,t(i))).eq.t(i)) then - call setei2(neigh(3,t(i)),t(i)) - else - call setei(t(i)) - endif - 10 continue -c -c perform local relaxations -c - call locrlx(t1) -c - return - end -c -c -------- divtri -c - subroutine divtri(t) - include 'commons' - integer t,t1,t2,stkpnt -c -c divide the triangle t by bisection -c -c construct the chain of triangles to be divided for compatability -c - stkpnt = 1 - stack(1) = t - t1 = t - 1 if (neigh(3,t1).le.0) go to 2 - if (neigh(3,neigh(3,t1)).eq.t1) go to 2 - stkpnt=stkpnt+1 - t1 = neigh(3,t1) - stack(stkpnt) = t1 - go to 1 - 2 continue -c -c debug output -c - if (outlev.ge.4) then - write(ioutpt,100) - write(ioutpt,101) (stack(i),i=1,stkpnt) - 100 format(' stack of triangles to divide is') - 101 format(13i6) - endif -c -c divide triangles in the chain -c - 3 if (stkpnt .eq. 0 .or. ierr .ne. 0 .or. - 1 (nvert .ge. ntarg .and. .not. unifrm)) go to 4 - t1 = stack(stkpnt) -c -c get the neighbor across from the peak of t1 -c - t2 = neigh(3,t1) -c -c divide t1 and t2 -c make sure the one with the error indicator comes first -c - if (coefei(iorder*iorder+1,1,t1).ne.0.) then - call dtpair(t1,t2) - else - call dtpair(t2,t1) - endif -c - stkpnt = stkpnt - 1 - go to 3 - 4 continue -c - return - end -c -c -------- adjds -c - subroutine adjds(t1,t2) - include 'commons' - integer t1,t2,level,i,nbr,v1,v2,v3,v4,i3,i4 -c -c adjust the data structures for vertices and triangles to -c divide the pair of triangles t1 and t2 -c -c debug output -c - if (outlev.ge.4) write(ioutpt,101) - 101 format(' adjust data structures for', - 1 ' vertices and triangles') -c -c set vertex numbers -c - v1 = vertex(1,t1) - v2 = vertex(2,t1) - v3 = vertex(3,t1) - v4 = vertex(3,t2) -c -c vertex data structures -c - nvert = nvert+1 - level = iabs(vrtlev(v3)) + 1 - if (level .gt. nlev) nlev = level - lvert(level) = lvert(level)+1 - vrtlev(nvert) = level - nextvt(nvert) = frstvt(level) - frstvt(level) = nvert - xvert(nvert) = (xvert(v1)+xvert(v2))/2. - yvert(nvert) = (yvert(v1)+yvert(v2))/2. -c -c define new nodes -c - call defnod(t1,t2) - do 9 i=nwndvt(nvert),nnode - inuse(i)=.false. - 9 continue -c -c triangle data structures -c - ntri = ntri + 2 -c - vertex(1,ntri-1) = v2 - vertex(2,ntri-1) = v3 - vertex(3,ntri-1) = nvert - vertex(1,ntri ) = v2 - vertex(2,ntri ) = v4 - vertex(3,ntri ) = nvert - vertex(2,t1 ) = v3 - vertex(3,t1 ) = nvert - vertex(2,t2 ) = v4 - vertex(3,t2 ) = nvert -c - neigh(1,ntri-1) = t1 - neigh(2,ntri-1) = ntri - neigh(3,ntri-1) = neigh(1,t1) - neigh(1,ntri ) = t2 - neigh(2,ntri ) = ntri-1 - neigh(3,ntri ) = neigh(1,t2) - neigh(3,t1 ) = neigh(2,t1) - neigh(2,t1 ) = t2 - neigh(1,t1 ) = ntri-1 - neigh(3,t2 ) = neigh(2,t2) - neigh(2,t2 ) = t1 - neigh(1,t2 ) = ntri -c -c correct neighbor of neighbors -c - nbr = neigh(3,ntri-1) - if (nbr .gt. 0) then - do 10 i=1,3 - if (neigh(i,nbr) .eq. t1) neigh(i,nbr) = ntri-1 - 10 continue - endif -c - nbr = neigh(3,ntri) - if (nbr .gt. 0) then - do 20 i=1,3 - if (neigh(i,nbr) .eq. t2) neigh(i,nbr) = ntri - 20 continue - endif -c -c set triangles next to vertices -c - do 30 i=8,1,-1 - if (tringl(i,v2).eq.t1) tringl(i,v2)=ntri-1 - if (tringl(i,v2).eq.t2) tringl(i,v2)=ntri - if (tringl(i,v3).eq.0) i3=i - if (tringl(i,v4).eq.0) i4=i - tringl(i,nvert)=0 - 30 continue - tringl(i3,v3) = ntri-1 - tringl(i4,v4) = ntri - tringl(1,nvert) = t1 - tringl(2,nvert) = t2 - tringl(3,nvert) = ntri-1 - tringl(4,nvert) = ntri -c -c remove t1 from error indicator lists -c - call elstrm(t1) -c -c debug output -c - if(outlev.ge.5) then - write(ioutpt,201) - write(ioutpt,300) t1,(vertex(i,t1),i=1,3), - 1 (neigh(i,t1),i=1,3) - write(ioutpt,300) t2,(vertex(i,t2),i=1,3), - 1 (neigh(i,t2),i=1,3) - write(ioutpt,300) ntri-1,(vertex(i,ntri-1),i=1,3), - 1 (neigh(i,ntri-1),i=1,3) - write(ioutpt,300) ntri,(vertex(i,ntri),i=1,3), - 1 (neigh(i,ntri),i=1,3) - write(ioutpt,301) nvert,vrtlev(nvert),xvert(nvert), - 1 yvert(nvert) - write(ioutpt,302) t1,(node(i,t1),i=1,nnodtr) - write(ioutpt,302) t2,(node(i,t2),i=1,nnodtr) - write(ioutpt,302) ntri-1,(node(i,ntri-1),i=1,nnodtr) - write(ioutpt,302) ntri,(node(i,ntri),i=1,nnodtr) - endif - 201 format(' new triangles, vertices, and neighbors are:') - 300 format(1x,i5,5x,3i5,5x,3i5) - 301 format(' new vertex ',i5,' has level',i4/ - 1 ' coordinates',2(2x,1pe15.8)) - 302 format(' nodes for triangle ',i5/10(10i6/)) -c - return - end -c -c -------- adjdsb -c - subroutine adjdsb(t) - include 'commons' - integer t,i,level,nbr,v1,v2,v3,i3 -c -c adjust the data structures for vertices and triangles to -c divide boundary triangle t -c - if (outlev .ge. 4) write(ioutpt,101) - 101 format(' adjust data structures for', - 1 ' vertices and triangles') -c -c set vertex numbers -c - v1 = vertex(1,t) - v2 = vertex(2,t) - v3 = vertex(3,t) -c -c vertex data structures -c - nvert = nvert + 1 - level = iabs(vrtlev(v3)) + 1 - if (level .gt. nlev) nlev=level - lbvert(level) = lbvert(level)+1 -c negate level to indicate dirichlet b.c. - vrtlev(nvert) = -level - nextvt(nvert) = frstvt(level) - frstvt(level) = nvert - xvert(nvert) = (xvert(v1) + xvert(v2))/2. - yvert(nvert) = (yvert(v1) + yvert(v2))/2. -c -c define new nodes -c - call defnod(t,0) - do 9 i=nwndvt(nvert),nnode - inuse(i)=.false. - 9 continue -c -c triangle data structures -c - ntri = ntri + 1 -c - vertex(1,ntri) = v2 - vertex(2,ntri) = v3 - vertex(3,ntri) = nvert - vertex(2,t ) = v3 - vertex(3,t ) = nvert -c - neigh(1,ntri) = t - neigh(2,ntri) = neigh(3,t) - neigh(3,ntri) = neigh(1,t) - neigh(3,t ) = neigh(2,t) - neigh(2,t ) = neigh(2,ntri) - neigh(1,t ) = ntri -c -c correct neighbor of neighbor -c - nbr = neigh(3,ntri) - if (nbr .gt. 0) then - do 10 i=1,3 - if (neigh(i,nbr) .eq. t) neigh(i,nbr) = ntri - 10 continue - endif -c -c set triangles next to vertices -c - do 30 i=8,1,-1 - if (tringl(i,v2).eq.t) tringl(i,v2)=ntri - if (tringl(i,v3).eq.0) i3=i - tringl(i,nvert)=0 - 30 continue - tringl(i3,v3) = ntri - tringl(1,nvert) = t - tringl(2,nvert) = ntri -c -c remove t from error indicator lists -c - call elstrm(t) -c -c debug output -c - if (outlev .ge. 5) then - write(ioutpt,201) - write(ioutpt,300) t,(vertex(i,t),i=1,3),(neigh(i,t),i=1,3) - write(ioutpt,300) ntri,(vertex(i,ntri),i=1,3), - 1 (neigh(i,ntri),i=1,3) - write(ioutpt,301) nvert,vrtlev(nvert),xvert(nvert), - 1 yvert(nvert) - write(ioutpt,302) t,(node(i,t),i=1,nnodtr) - write(ioutpt,302) ntri,(node(i,ntri),i=1,nnodtr) - endif - 201 format(' new triangles, vertices, and neighbors are:') - 300 format(1x,i5,5x,3i5,5x,3i5) - 301 format(' new vertex ',i5,' has level',i4/' coordinates', - 1 2(2x,1pe15.8)) - 302 format(' nodes and bndnod for triangle ',i5/10(10i6/)) -c - return - end -c -c -------- defnod -c - subroutine defnod(t1,t2) - include 'commons' - integer i,j,t1,t2,limit,bctype - real bcrhs,cu -c -c define the nodes for the descendents of the triangle pair t1,t2 -c -c set first new node associated with the new vertex -c - nwndvt(nvert)=nnode+1 -c -c set old nodes associated with the new vertex -c - do 10 i=1,nnodtr - olndvt(i,nvert)=node(i,t1) - 10 continue - if (t2.ne.0) then - do 20 i=iorder+1,nnodtr - olndvt(nnodtr-iorder+i,nvert)=node(i,t2) - 20 continue - endif -c -c default new nodes as interior, fix boundary later -c - if (t2.ne.0) then - limit=nnodev - else - limit=nnodvb - endif - do 30 i=1,limit - bndnod(nnode+i)=.false. - 30 continue -c -c descendents of t1 -c - do 40 i=1,nnodtr - hldnod(i) = node(i,t1) - 40 continue -c -c nodes for the new t1 -c - do 50 i=1,nnodtr - j=renum(1,i) - if (j.gt.0) then - node(i,t1)=hldnod(j) - else - node(i,t1)=nnode-j - endif - 50 continue -c -c nodes for the other part of t1 -c - do 60 i=1,nnodtr - j=renum(2,i) - if (j.gt.0) then - node(i,ntri+1)=hldnod(j) - else - node(i,ntri+1)=nnode-j - endif - 60 continue -c -c descendents of t2. If t2=0 then t1 is boundary -c - if (t2 .ne. 0) then -c - do 70 i=1,nnodtr - hldnod(i) = node(i,t2) - 70 continue -c -c nodes for the new t2 -c - do 80 i=1,nnodtr - j=renum(3,i) - if (j.gt.0) then - node(i,t2)=hldnod(j) - else - node(i,t2)=nnode-j - endif - 80 continue -c -c nodes for the other part of t2 -c - do 90 i=1,nnodtr - j=renum(4,i) - if (j.gt.0) then - node(i,ntri+2)=hldnod(j) - else - node(i,ntri+2)=nnode-j - endif - 90 continue - nnode=nnode+nnodev -c - else -c -c if t2 is 0, mark boundary nodes and copy boundary condition -c from error indicator problem -c -c@@@ call is just to get bctype - call bcond(xvert(vertex(1,t1)),yvert(vertex(1,t1)), - . -neigh(3,t1),cu,bcrhs,bctype) - if (bctype.eq.1) then - do 110 i=1,nnodvb - if (nbasch(i).eq.iorder) then - bndnod(nnode+i)=.true. - u(nnode+i)=bcei(i,t1) - endif - 110 continue - endif - nnode=nnode+nnodvb -c - endif -c - return - end -c -c -------- oldip -c - subroutine oldip(t1,t2) - include 'commons' - integer t1,t2,limit,i,j,t -c -c compute the inner products of the old nodes over the 2 -c old triangles t1 and t2 and subtract from the matrix and rs. -c if t2<=0 (boundary) 1 triangle -c -c debug output -c - if (outlev.ge.4) write(ioutpt,101) - 101 format(' subtract old inner products') -c - if (t2.le.0) then - limit=1 - else - limit=2 - endif -c -c for each of the triangles -c - do 20 i=1,limit -c -c pick triangle -c - if (i.eq.1) then - t=t1 - else - t=t2 - endif -c -c do quadratures -c - call quad(t,1) -c -c update right side -c - do 10 j=1,naddrs - rs(rowrs(j))=rs(rowrs(j))-addrs(j) - if (outlev.ge.5) write(ioutpt,102) rowrs(j),rs(rowrs(j)) - 10 continue -c -c negate add to subtract -c - do 15 j=1,nadd - add(j)=-add(j) - 15 continue -c -c update matrix -c - call addmat -c - 20 continue - 102 format(' rowrs(',i6,') <-- ',1pe15.8) -c - return - end -c -c -------- newip -c - subroutine newip(t1,t2) - include 'commons' - integer t1,t2,limit,i,j,t -c -c compute the new inner products of the old nodes over the 4 -c new triangles (t1,t2,ntri-1 and ntri) and add to matrix and rs -c if t2<=0 (boundary) 2 triangles (t1,ntri) -c -c debug output -c - if (outlev.ge.4) write(ioutpt,101) - 101 format(' add new inner products') -c - if (t2.le.0) then - limit=2 - else - limit=4 - endif -c -c for each of the triangles -c - do 20 i=1,limit -c -c pick triangle -c - if (i.eq.1) then - t=t1 - elseif (i.eq.2) then - t=ntri - elseif (i.eq.3) then - t=t2 - else - t=ntri-1 - endif -c -c do quadratures -c - call quad(t,2) -c -c update right side -c - do 10 j=1,naddrs - rs(rowrs(j))=rs(rowrs(j))+addrs(j) - if (outlev.ge.5) write(ioutpt,102) rowrs(j),rs(rowrs(j)) - 10 continue -c -c update matrix -c - call addmat -c - 20 continue - 102 format(' rowrs(',i6,') <-- ',1pe15.8) -c - return - end - -c -c -------- neweq -c - subroutine neweq(t1,t2) - include 'commons' - integer t1,t2,i,j,node1 -c -c set new rows of matrix (from dividing triangle pair t1,t2) -c -c debug output -c - if (outlev.ge.4) write(ioutpt,101) t1 - 101 format(' set new equations from dividing triangle ',i6) -c -c set idcoef for new rows -c - call setidn(t1) -c -c change idcoef for effected old rows -c - call convid(nvert,1) -c -c copy rs and coef from error indicator problem -c - node1=nwndvt(nvert) - do 30 i=node1,nnode - rs(i)=rsei(i-node1+1,t1) - do 10 j=1,iorder*iorder - coef(j,i)=coefei(j,i-node1+1,t1) - 10 continue - do 20 j=iorder*iorder+1,mxidlo - coef(j,i)=0. - 20 continue - coef(mxidlo+1,i)=coefei(iorder*iorder+1,i-node1+1,t1) - 30 continue -c -c get values for coef for new node - new node inner products -c from bloke2 -c - do 50 i=node1,nnode - do 50 j=1,mxidlo - if (idcoef(j,i).ge.node1) then - coef(j,i)=bloke2(idcoef(j,i)-node1+1,i-node1+1,t1) - endif - 50 continue -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,102) - do 110 i=nwndvt(nvert),nnode - write(ioutpt,103) (idcoef(j,i),coef(j,i),j=1,mxidlo), - 1 i,coef(mxidlo+1,i) - write(ioutpt,104) rs(i) - write(ioutpt,105) (idcoef(j,i),j=mxidlo+1,mxidup) - 110 continue - endif - 102 format(' new equations are ') - 103 format(3(i6,1pe15.8)) - 104 format(' rs ',1pe15.8) - 105 format(11i6) -c -c local relaxation at new nodes -c - call rlxred(nvert) -c -c debug output -c - if (outlev.ge.5) write(ioutpt,100) (i,u(i),i=nwndvt(nvert),nnode) - 100 format(' initial guess for node ',i6,2x,1pe15.8) -c - return - end -c -c -------- setidn -c - subroutine setidn(t1) - include 'commons' - integer node1,i,j,limit,t,tri,nod1,nod2,c,t1 -c -c set idcoef for new rows -c -c copy idcoef for old node - new node inner products from -c error indicator problem -c - node1=nwndvt(nvert)-1 - if (vrtlev(nvert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 30 i=1,limit - do 10 j=1,iorder*iorder - idcoef(j,node1+i)=idcoei(j,i,t1) - 10 continue - do 20 j=iorder*iorder+1,mxidup - idcoef(j,node1+i)=0 - 20 continue - 30 continue -c -c set idcoef for new node -- new node inner products -c - if (iorder.gt.2) then - if (vrtlev(nvert).gt.0) then - limit=4 - else - limit=2 - endif -c - do 130 t=1,limit - tri=tringl(t,nvert) - do 120 i=1,nnodtr - nod1=node(i,tri) - if (nod1.gt.node1) then - do 110 j=1,nnodtr - nod2=node(j,tri) - if (nod2.gt.node1 .and. nod2.ne.nod1) then - if(nod2.lt.nod1) then - c=1 - 1 if (idcoef(c,nod1).eq.0 .or. - 1 idcoef(c,nod1).eq.nod2) go to 2 - if (c.ge.mxidlo) then - write(ioutpt,101) - stop - 101 format(' ********FATAL ERROR********'// - 1 ' ran out of room in idcoef', - 1 ' lower part in setidn') - endif - c=c+1 - go to 1 - 2 continue - idcoef(c,nod1)=nod2 - else - c=mxidlo+1 - 3 if (idcoef(c,nod1).eq.0 .or. - 1 idcoef(c,nod1).eq.nod2) go to 4 - if (c.ge.mxidup) then - write(ioutpt,102) - stop - 102 format(' ********FATAL ERROR********'// - 1 ' ran out of room in idcoef', - 1 ' upper part in setidn') - endif - c=c+1 - go to 3 - 4 continue - idcoef(c,nod1)=nod2 - endif - endif - 110 continue - endif - 120 continue - 130 continue - endif -c - return - end -c -c -------- addmat -c - subroutine addmat - include 'commons' - integer i,j,c,r -c -c add nadd values in add to the matrix in -c positions given by (row,col) -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) nadd - endif -c -c for each change, find the idcoef with col(i) in row(i) -c and add add(i) -c -c begin the search for the column of coef with the position -c that the previous value was in, and wrap around -c - lastc=1 - do 20 i=1,nadd - r=row(i) - if (row(i).eq.col(i)) then -c special case for diagonal entry - c=mxidlo+1 - lastc=1 - else - c=0 -c search for column - j=lastc - jfirst=j-1 - if (jfirst.le.0) jfirst=mxidlo - 1 if (j.eq.jfirst .or. idcoef(j,r).eq.col(i)) go to 2 - if (j.ge.mxidlo .or. idcoef(j,r).eq.0) j=0 - j=j+1 - go to 1 - 2 continue -c verify that the column was found - if (idcoef(j,r).eq.col(i)) c=j - lastc=c - endif - if(c.eq.0) go to 999 -c add - coef(c,r) = coef(c,r) + add(i) -c - 20 continue - return -c -c error that should not occur - couldnt find location -c - 999 write(ioutpt,100) col(i),row(i) - stop -c - 100 format(' ********FATAL ERROR********'// - 1 ' couldnt find column ',i10,' in row ',i10, - 2 ' while adding a value to the matrix') - 200 format(' add ',i3,' values to matrix') -c - end -c -c -------- quad -c - subroutine quad(t,itype) - include 'commons' - integer t,itype,i,j,k,sub,qp,node1,node2,bctype,nlist(ndord) - real x1,x2,x3,y1,y2,y3,det,area,dz1dx, - 1 dz2dx,dz3dx,dz1dy,dz2dy,dz3dy,qpx,qpy, - 2 db1dx,db2dx,db1dy,db2dy - real bcrhs,cu - real qpf,qpp,qpq,qpr - logical bndsid(3) -c -c compute integrals over triangle t -c -c itype determines which integrals to compute -c = 1 all -c = 2 old node - right side and old node - old node -c = 3 new node - right side, new node - new node and -c new node - old node -c -c the nadd values for the matrix in add go in positions (row,col) -c the naddrs values for the right side in addrs go rows rowrs -c -c compute linear transformation from reference triangle -c - x1=xvert(vertex(1,t)) - x2=xvert(vertex(2,t)) - x3=xvert(vertex(3,t)) - y1=yvert(vertex(1,t)) - y2=yvert(vertex(2,t)) - y3=yvert(vertex(3,t)) -c - det = x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2) - area = abs(det/2.) -c - dz1dx = (y2-y3)/det - dz2dx = (y3-y1)/det - dz3dx = (y1-y2)/det - dz1dy = (x3-x2)/det - dz2dy = (x1-x3)/det - dz3dy = (x2-x1)/det -c -c determine which sides of the triangle are boundary with -c non-Dirichlet boundary conditions -c - do 5 i=1,3 - bndsid(i) = .false. - if (neigh(i,t).le.0) then -c@@@ call is just to get bctype - call bcond(x1,y1,-neigh(i,t),cu,bcrhs,bctype) - if (bctype.ne.1) bndsid(i) = .true. - endif - 5 continue -c -c in local numbering, determine which bases to use and -c initialize right side at 0 -c - if (itype.eq.1) then -c compute all inner products - do 10 i=1,nnodtr - rowrs(i)=i - addrs(i)=0. - 10 continue - naddrs=nnodtr -c - elseif (itype.eq.2) then -c -c only inner products of old nodes with old nodes - naddrs=0 - sub=0 - do 30 i=1,(iorder+1)/2 - do 20 j=1,iorder-2*i+2 - sub=sub+1 - naddrs=naddrs+1 - rowrs(naddrs)=sub - addrs(naddrs)=0. - 20 continue - sub=sub+iorder-2*i+1 - 30 continue -c - else -c -c only new nodes with new nodes and old nodes - naddrs=0 - sub=iorder - do 50 i=1,iorder/2 - do 40 j=1,iorder-2*i+1 - sub=sub+1 - naddrs=naddrs+1 - rowrs(naddrs)=sub - addrs(naddrs)=0. - 40 continue - sub=sub+iorder-2*i - 50 continue -c - endif -c -c set row and col for all combinations of nodes used and -c initalize inner products at 0 -c - nadd=0 - do 60 i=1,naddrs - do 60 j=i,naddrs - nadd=nadd+1 - row(nadd)=rowrs(j) - col(nadd)=rowrs(i) - add(nadd)=0. - 60 continue -c -c for itype=3, extend row and col to include old node-new node -c - if (itype.eq.3) then - sub=0 - do 90 i=1,(iorder+1)/2 - do 80 j=1,iorder-2*i+2 - sub=sub+1 - do 70 k=1,naddrs - nadd=nadd+1 - row(nadd)=rowrs(k) - col(nadd)=sub - add(nadd)=0. - 70 continue - 80 continue - sub=sub+iorder-2*i+1 - 90 continue - endif -c -c -c compute integrals -c -c -c for each quadrature point . . . -c - do 130 qp=1,nqpt -c -c compute x-y coordinate of quadrature point -c - qpx = x1*quadpt(1,qp) + x2*quadpt(2,qp) + x3*quadpt(3,qp) - qpy = y1*quadpt(1,qp) + y2*quadpt(2,qp) + y3*quadpt(3,qp) -c -c evaluate pde coefficents and right side at quadrature point -c - call pde(qpx,qpy,qpp,qpq,qpr,qpf) -c -c for each basis pair . . . -c - do 110 i=1,nadd -c -c compute derivatives of bases wrt x and y at quadrature point -c - db1dx = qpdbdz(1,row(i),qp)*dz1dx - 1 + qpdbdz(2,row(i),qp)*dz2dx - 2 + qpdbdz(3,row(i),qp)*dz3dx - db1dy = qpdbdz(1,row(i),qp)*dz1dy - 1 + qpdbdz(2,row(i),qp)*dz2dy - 2 + qpdbdz(3,row(i),qp)*dz3dy - db2dx = qpdbdz(1,col(i),qp)*dz1dx - 1 + qpdbdz(2,col(i),qp)*dz2dx - 2 + qpdbdz(3,col(i),qp)*dz3dx - db2dy = qpdbdz(1,col(i),qp)*dz1dy - 1 + qpdbdz(2,col(i),qp)*dz2dy - 2 + qpdbdz(3,col(i),qp)*dz3dy -c -c contribution of this quadrature point to this integral -c - add(i) = add(i) + - 1 quadw(qp) * ( qpp*db1dx*db2dx + qpq*db1dy*db2dy - 2 + qpr*qpbas(row(i),qp)*qpbas(col(i),qp) ) -c - 110 continue -c -c right side integral -c - do 120 i=1,naddrs - addrs(i) = addrs(i) + - 1 quadw(qp) * qpf * qpbas(rowrs(i),qp) - 120 continue - 130 continue -c -c finish integrals by multiplying by the area of the triangle -c - do 140 i=1,naddrs - addrs(i) = addrs(i)*area - 140 continue - do 150 i=1,nadd - add(i) = add(i)*area - 150 continue -c -c compute boundary integrals -c - do 155 i=1,3 - if (bndsid(i)) then -c make list of nodes on this boundary side - if (i.eq.1) then - do 151 j=1,iorder - nlist(j) = renum(2,j) - 151 continue - xx1 = x2 - yy1 = y2 - xx2 = x3 - yy2 = y3 - elseif (i.eq.2) then - do 152 j=1,iorder - nlist(j) = renum(1,j) - 152 continue - xx1 = x1 - yy1 = y1 - xx2 = x3 - yy2 = y3 - else - do 153 j=1,iorder - nlist(j) = j - 153 continue - xx1 = x1 - yy1 = y1 - xx2 = x2 - yy2 = y2 - endif -c - call quadb(nlist,xx1,yy1,xx2,yy2,t,i) - endif - 155 continue -c -c change local numbering of nodes to global numbering -c - if (itype.ne.3) then -c - do 160 i=1,naddrs - rowrs(i)=node(rowrs(i),t) - 160 continue -c - do 170 i=1,nadd - node1=node(row(i),t) - node2=node(col(i),t) - if (node1.gt.node2) then - row(i)=node1 - col(i)=node2 - else - row(i)=node2 - col(i)=node1 - endif - 170 continue -c - endif -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) t - 200 format(' integrals over triangle ',i6/' row,col,matrix') - write(ioutpt,300) (row(i),col(i),add(i),i=1,nadd) - 300 format(1x,2i6,1pe15.8) - write(ioutpt,400) - 400 format(' row,right side') - write(ioutpt,500) (rowrs(i),addrs(i),i=1,naddrs) - 500 format(1x,i6,1pe15.8) - endif -c - return - end -c -c -------- quadb -c - subroutine quadb(nlist,x1,y1,x2,y2,t,iside) - include 'commons' - integer nlist(*),t - real x1,y1,x2,y2 - integer incl1(ndord6),incl2(ndord6),inclrs(ndord3) - integer i,j,qp,bctype - real qpx,qpy,slen - real cu,bcrhs -c -c compute boundary integrals -c -c x1,y1,x2,y2 are the vertices of the triangle on the boundary side -c t is the triangle -c nlist is the local numbering of the nodes on the boundary side -c - slen = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) -c -c determine which of the nadd matrix changes involve the boundary integral -c - do 30 j=1,nadd - incl1(j) = 0 - do 10 i=1,iorder - if (row(j).eq.nlist(i)) incl1(j) = i - 10 continue - incl2(j) = 0 - do 20 i=1,iorder - if (col(j).eq.nlist(i)) incl2(j) = i - 20 continue - 30 continue -c -c same for right side -c - do 50 j=1,naddrs - inclrs(j) = 0 - do 40 i=1,iorder - if (rowrs(j).eq.nlist(i)) inclrs(j) = i - 40 continue - 50 continue -c -c for each quadrature point -c - do 80 qp=1,nqptb -c -c set x,y coordinate of quadrature point -c - qpx = x1*qptb(qp) + x2*(1.-qptb(qp)) - qpy = y1*qptb(qp) + y2*(1.-qptb(qp)) -c -c evaluate boundary condition -c - call bcond(qpx,qpy,-neigh(iside,t),cu,bcrhs,bctype) -c -c contributions from this quadrature point -c - do 60 i=1,nadd - if (incl1(i).ne.0 .and. incl2(i).ne.0) then - add(i) = add(i) + qwtb(qp)* - 1 (cu*qpbasb(incl1(i),qp)*qpbasb(incl2(i),qp))*slen - endif - 60 continue -c - do 70 i=1,naddrs - if (inclrs(i).ne.0) then - addrs(i) = addrs(i) - 1 + qwtb(qp)*bcrhs*qpbasb(inclrs(i),qp)*slen - endif - 70 continue - 80 continue - return - end -c -c -------- solve -c - subroutine solve - include 'commons' - real oemaxv,oemaxn,oemaxq,oeenrg - external clock -c -c solve the system of equations coef * u = rs by ncyc multigrid V-cycles -c u is the coefficient vector for the nodal basis -c u is used as initial guess -c -c debug output -c - if (outlev.ge.2) write(ioutpt,100) - 100 format(/' begin solution ') -c -c start timing -c - timesl = clock(0) -c -c initial errors -c - if (outlev.ge.3) then - call errors(oemaxv,oemaxn,oemaxq,oeenrg) - write(ioutpt,101) oemaxv,oemaxn,oemaxq,oeenrg - endif - 101 format(/' errors in initial guess'/ - 1 ' max norm at vertices ',1pe15.8/ - 2 ' max norm at nodes ',1pe15.8/ - 3 ' max norm at quad pts ',1pe15.8/ - 4 ' continuous energy norm ',1pe15.8/) -c -c solve approximatly using ncyc cycles -c - do 902 loop=1,ncyc -c - if (outlev.ge.3) write(ioutpt,201) - 201 format(' begin multigrid cycle') -c -c v-cycle -c - if (nlev.ne.1) then - if (outlev.ge.3) write(ioutpt,301) - 301 format(' restrictions') - do 10 l=nlev,2,-1 - call relax(l,nu1) - call restrc(l) - 10 continue - endif -c - call exsolv -c - if (nlev.ne.1) then - if (outlev.ge.3) write(ioutpt,303) - 303 format(' prolongations') - do 20 l=2,nlev - call prolon(l) - call relax(l,nu2) - 20 continue - endif -c - if (outlev.ge.3) write(ioutpt,202) - 202 format(' multigrid cycle complete') - if (outlev.ge.3) call errred(loop,oemaxv,oemaxn, - 1 oemaxq,oeenrg) - 902 continue -c -c finish timing -c - timesl = clock(0) - timesl - timest = timest + timesl -c -c write data files for run time graphics, check the menu, and -c update the plots -c - if (outlev.ge.2 .and. grafic) write(ioutpt,200) - 200 format(' begin update to solution graphics') -c - if (grafic) then - call filtri - call filsol - if (menuon) call mn2mg - call grfsol - endif - if (outlev.ge.2 .and. grafic) write(ioutpt,300) - 300 format(' graphics complete') -c - call outsol -c - return - end -c -c -------- restrc -c - subroutine restrc(lev) - include 'commons' - integer lev,vert,limit,node1,i,nod,k -c -c compute residual on grid lev and restrict to grid lev-1 -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) lev,lev-1 - 100 format(' restrict from level',i4,' to level',i4) -c -c convert system from nodal basis to 2-level hierarchical -c - call stmats(lev,-1) - call svec(u,lev,-1) - call stvec(rs,lev,-1) -c -c compute part of residual due to level lev -c -c debug output -c - if (outlev .ge. 4) write(ioutpt,200) lev - 200 format(' compute residual due to level',i4) -c - vert=frstvt(lev) - 1 if (vert .eq. -1) go to 2 - node1=nwndvt(vert) - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 10 i=1,limit - nod=node1-1+i - k=1 - 11 if (idcoef(k,nod).eq.0 .or. idcoef(k,nod).ge.node1) - 1 go to 12 - rs(idcoef(k,nod)) = rs(idcoef(k,nod)) - 1 - coef(k,nod)*u(nod) - k=k+1 - go to 11 - 12 continue - 10 continue - vert = nextvt(vert) - go to 1 - 2 continue -c - return - end -c -c -------- prolon -c - subroutine prolon(lev) - include 'commons' - integer lev,vert,limit,node1,nod,k,i -c -c prolongate from grid lev-1 to grid lev -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) lev-1,lev - 100 format(' prolongate from level',i4,' to level',i4) - if (outlev.ge.4) write(ioutpt,101) lev - 101 format(' compute residual due to level ',i6) -c -c compute part of residual due to level lev -c and subtract it off -c - vert=frstvt(lev) - 1 if (vert .eq. -1) go to 2 - node1=nwndvt(vert) - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 10 i=1,limit - nod=node1-1+i - k=1 - 11 if (idcoef(k,nod).eq.0 .or. idcoef(k,nod).ge.node1) - 1 go to 12 - rs(idcoef(k,nod)) = rs(idcoef(k,nod)) - 1 + coef(k,nod)*u(nod) - k=k+1 - go to 11 - 12 continue - 10 continue - vert = nextvt(vert) - go to 1 - 2 continue -c -c -c convert system from 2-level hierarchical basis to nodal basis -c - call stmats(lev,1) - call svec(u,lev,1) - call stvec(rs,lev,1) -c - return - end -c -c -------- svec -c - subroutine svec(vec,lev,dir) - include 'commons' - real vec(*) - integer lev,vert,node1,nodei,i,j,oldnod,dir,limit -c -c convert vector vec from a 2-level coefficent vector to a -c nodal coefficient vector on level lev (dir=1) -c or nodal to 2-level (dir=-1) -c i.e. multiply vec by s (or s inverse) where s is the matrix that -c converts a 2-level basis to a nodal basis on level lev -c - if (lev.le.1) return -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) dir,lev - 100 format(' change basis of solution vector in direction',i3, - 1 ' on level',i4) -c - vert=frstvt(lev) - 1 if (vert.eq.-1) go to 2 - node1=nwndvt(vert)-1 - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 20 i=1,limit - nodei=node1+i - do 10 j=1,nbasch(i) - oldnod=olndvt(ibasch(j,i),vert) - vec(nodei)=vec(nodei)+dir*cbasch(j,i)*vec(oldnod) - 10 continue - 20 continue - vert = nextvt(vert) - go to 1 - 2 continue - 30 continue -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) - call prvec(vec,1,lev) - 200 format(' converted vector') - endif -c - return - end -c -c -------- stvec -c - subroutine stvec(vec,lev,dir) - include 'commons' - real vec(*) - integer vert,lev,node1,nodei,i,j,oldnod,dir,limit -c -c change a right side vector, vec, for converting a nodal basis -c on level lev to a 2-level hierarchical basis (dir=-1) -c or 2-level to nodal (dir=1) -c i.e. multiply the vector vec by s transpose (or s inverse transpose) -c where s is the matrix that converts a 2-level basis -c to a nodal basis on level lev -c - if (lev.le.1) return -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) dir,lev - 100 format(' change basis of right side vector in direction',i3, - 1 ' on level ',i4) -c - vert=frstvt(lev) - 1 if (vert .eq. -1) go to 2 -c - node1=nwndvt(vert)-1 - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 20 i=1,limit - nodei=node1+i - do 10 j=1,nbasch(i) - oldnod=olndvt(ibasch(j,i),vert) - vec(oldnod)=vec(oldnod)-dir*cbasch(j,i)*vec(nodei) - 10 continue - 20 continue - vert = nextvt(vert) - go to 1 - 2 continue -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) - call prvec(vec,1,lev) - 200 format(' converted vector') - endif -c - return - end -c -c -------- stmats -c - subroutine stmats(lev,dir) - include 'commons' - integer lev,dir,i,j,k,vert,limit,node1,nodei, - 1 oldnod,sub,nodek - real c -c -c change the matrix from nodal basis to 2-level hierarchical -c basis (dir=-1) on level lev, or 2-level to nodal (dir=1) -c i.e. compute s transpose * matrix * s where s is the matrix -c that converts a 2-level hierarchical basis to a nodal basis -c on level lev (or s inverse transpose * matrix * s inverse) -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) lev,dir - 100 format(' change basis of matrix on level',i4, - 1 ' in direction ',i2) -c - vert=frstvt(lev) -c - 1 if (vert.eq.-1) go to 2 -c - node1=nwndvt(vert)-1 - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - do 20 i=1,limit - nodei=node1+i - do 10 j=1,nbasch(i) - oldnod=olndvt(ibasch(j,i),vert) - c=-dir*cbasch(j,i) - sub=1 -c - 3 if (idcoef(sub,nodei).eq.oldnod) go to 4 - sub=sub+1 - if (sub.gt.mxidlo) then - write(ioutpt,101) oldnod,nodei - stop - endif - go to 3 - 4 continue -c - add(1)=2.*c*coef(sub,nodei) - 1 +c*c*coef(mxidlo+1,nodei) - row(1)=oldnod - col(1)=oldnod - nadd=1 -c - k=1 - 5 if (k.gt.mxidlo.or.idcoef(k,nodei).eq.0) go to 6 - nodek=idcoef(k,nodei) - nadd=nadd+1 - add(nadd)=c*coef(k,nodei) - if (nodek.gt.oldnod) then - row(nadd)=nodek - col(nadd)=oldnod - elseif (nodek.lt.oldnod) then - row(nadd)=oldnod - col(nadd)=nodek - else - nadd=nadd-1 - endif - k=k+1 - go to 5 - 6 continue -c - nadd=nadd+1 - add(nadd)=c*coef(mxidlo+1,nodei) - row(nadd)=nodei - col(nadd)=oldnod -c - k=mxidlo+1 - 7 if (k.gt.mxidup.or.idcoef(k,nodei).eq.0) go to 8 - nodek=idcoef(k,nodei) - if (nodek .le. node1+limit) then - sub=1 - 11 if (idcoef(sub,nodek).eq.nodei) go to 12 - sub=sub+1 - if (sub.gt.mxidlo) then - write(ioutpt,101) nodei,nodek - stop - endif - go to 11 - 12 continue - nadd=nadd+1 - add(nadd)=c*coef(sub,nodek) - row(nadd)=nodek - col(nadd)=oldnod - endif - k=k+1 - go to 7 - 8 continue -c - call addmat - 10 continue - 20 continue - call convid(vert,dir) -c - vert=nextvt(vert) - go to 1 - 2 continue -c - 101 format(/' ********FATAL ERROR********'// - 1 ' couldnt find column ',i10,' in row ',i10, - 2 ' in subroutine stmats ') -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) - call outmat - 200 format(' matrix after basis change') - endif -c - return - end -c -c -------- convid -c - subroutine convid(vert,dir) - include 'commons' - integer limit,i,j,count,vert,dir,inc -c -c convert upper triangular part of idcoef form 2-level to nodal -c (dir=1) or nodal to 2-level (dir=-1) at vertex vert -c - if (vrtlev(vert).gt.0) then - limit=3 - else - limit=1 - endif -c -c for each of the two macrotriangles (1 if boundary) -c - do 99 i=1,limit,2 -c -c make a list of the old vertices in the half with the oldest -c vertex which are not in both halves -c - count=0 - do 10 j=1,nnodtr - if (renum(i,j).gt.0) then - count=count+1 - if (i.eq.1 .or. renum(i,j).le.iorder) then - lolnd1(count)=olndvt(renum(i,j),vert) - else - lolnd1(count)=olndvt(renum(i,j)+nnodtr-iorder,vert) - endif - elseif (renum(i,j-1).gt.0) then - count=count-1 - endif - 10 continue - if (renum(i,nnodtr).gt.0) count=count-1 -c -c make a list of the old vertices in the other half which are -c not in both halves -c - count=0 - do 20 j=1,nnodtr - if (renum(i+1,j).gt.0) then - count=count+1 - lolnd2(count)=olndvt(renum(i+1,j),vert) - if (i.eq.1 .or. renum(i+1,j).le.iorder) then - lolnd2(count)=olndvt(renum(i+1,j),vert) - else - lolnd2(count)=olndvt(renum(i+1,j)+nnodtr-iorder,vert) - endif - elseif (renum(i+1,j-1).gt.0) then - count=count-1 - endif - 20 continue - if (renum(i+1,nnodtr).gt.0) count=count-1 -c -c make a list of the new nodes in the half with the oldest vertex -c - count=0 - do 30 j=1,nnodtr - if (renum(i,j).lt.0) then - count=count+1 - lnewnd(count)=nwndvt(vert)-1-renum(i,j) - endif - 30 continue -c -c for each node in old node list 1, replace old node list 2 by -c new node list in idcoef (or other way depending on direction) -c - do 40 j=1,count - if (dir.eq.1) then - call replac(lolnd2,lnewnd,count,lolnd1(j)) - else - call replac(lnewnd,lolnd2,count,lolnd1(j)) - endif - 40 continue -c -c make a list of the new nodes in the half without the oldest vertex -c - count=0 - do 50 j=1,nnodtr - if (renum(i+1,j).lt.0) then - count=count+1 - lnewnd(count)=nwndvt(vert)-1-renum(i+1,j) - endif - 50 continue -c -c for each node in old node list 2, replace old node list 1 by -c new node list in idcoef (or other way depending on direction) -c - do 60 j=1,count - if (dir.eq.1) then - call replac(lolnd1,lnewnd,count,lolnd2(j)) - else - call replac(lnewnd,lolnd1,count,lolnd2(j)) - endif - 60 continue -c -c make a list of new nodes in both halves of the triangle -c - count=0 - do 70 j=1,nnodtr - if (renum(i,j).lt.0) then - count=count+1 - lnewnd(count)=nwndvt(vert)-1-renum(i,j) - endif - 70 continue - do 80 j=1,nnodtr - if(renum(i+1,j).lt.0) then - count=count+1 - lnewnd(count)=nwndvt(vert)-1-renum(i+1,j) - elseif (j.ne.1) then - if (renum(i+1,j-1).lt.0) count=count-1 - endif - 80 continue - if (renum(i+1,nnodtr).lt.0) count=count-1 -c -c make a list of 0's of the same length -c - do 90 j=1,count - lolnd1(j)=0 - 90 continue -c -c for each node on boundary between the 2 halves, replace 0's -c by new nodes (or other way depending on direction) -c - j=iorder - inc=1 - 91 if (i.eq.1 .or. renum(i,j).le.iorder) then - nod=olndvt(renum(i,j),vert) - else - nod=olndvt(renum(i,j)+nnodtr-iorder,vert) - endif - if (dir.eq.1) then - call replac(lolnd1,lnewnd,count,nod) - else - call replac(lnewnd,lolnd1,count,nod) - endif - if (j.eq.nnodtr .or. j.eq.nnodtr-1) go to 92 - j=j+2*(iorder-inc)-1 - inc=inc+2 - go to 91 - 92 continue -c - 99 continue -c - return - end -c -c -------- replac -c - subroutine replac(nodout,nodin,nrepl,irow) - include 'commons' - integer nodout(*),nodin(*),irow,nrepl,i,no,ni, - 1 icol,col2 -c -c replace the nrepl nodes in nodout in the upper triangular -c part of idcoef in row row by the nodes in nodin -c 0's are kept at the end of a row of idcoef -c nodes less than row are replaced by 0 -c - if (outlev.ge.5) then - write(ioutpt,200) irow - write(ioutpt,300) (nodout(i),i=1,nrepl) - write(ioutpt,400) (nodin(i),i=1,nrepl) - 200 format(' replace upper idcoef for row ',i6) - 300 format(' remove ',14i6) - 400 format(' add ',14i6) - endif -c - do 99 i=1,nrepl - no=nodout(i) - if (no.le.irow) no=0 - ni=nodin(i) - if (ni.le.irow) ni=0 -c -c if nodout=0, make sure nodin is not already there -c - if (no.eq.0) then - do 9 icol=mxidlo+1,mxidup - if (idcoef(icol,irow).eq.ni) ni=0 - 9 continue - endif - if(no.eq.0 .and. ni.eq.0) go to 99 -c -c find nodout -c - icol=mxidlo+1 - 1 if (idcoef(icol,irow).eq.no) go to 2 - icol=icol+1 -c if nodout is not found, the replacement was already done - if (icol.gt.mxidup) go to 99 - go to 1 - 2 continue -c -c replace with nodin -c - if (ni.ne.0) then - idcoef(icol,irow)=ni - else -c -c put 0 nodin at end -c - col2=icol - 3 if (idcoef(col2,irow).eq.0) go to 4 - col2=col2+1 - if (col2.le.mxidup) go to 3 - 4 continue - idcoef(icol,irow)=idcoef(col2-1,irow) - idcoef(col2-1,irow)=0 - endif -c - 99 continue -c - return - end -c -c -------- initxs -c - subroutine initxs - include 'commons' - integer i,j,id - real xnode(ndrow0),ynode(ndrow0) -c -c initialization for exact solve on coarsest grid by linpack band -c - if (outlev.ge.3) write(ioutpt,104) - 104 format(' begin initializations for coarse grid exact solve') -c -c verify enough space is allocated for coarse grid matrix rows -c - if (nnode0 .gt. ndrow0) then - write(ioutpt,101) ndrow0,nnode0 - 101 format(/'FATAL ERROR -- not enough room for coarsest grid', - 1 ' factors for exact solve'/ - 2 'Current allocation is ndrow0 = ',i6/ - 3 'Requires ndrow0 >= ',i6/ - 4 'Change parameter statement in commons') - stop - endif -c - if (outlev.ge.4) write(ioutpt,105) - 105 format(' determine reordering of coarse grid nodes') -c -c determine reordering of coarse grid nodes to reduce bandwidth -c -c initially the order is the order in the grid; also set flag in xnode -c - do 10 i=1,nnode0 - ipvt1(i)=i - xnode(i)=-999. - 10 continue -c -c compute node coordinates -c - delta = 1./float(iorder-1) - do 40 itri=1,ntri - x1 = xvert(vertex(1,itri)) - x2 = xvert(vertex(2,itri)) - x3 = xvert(vertex(3,itri)) - y1 = yvert(vertex(1,itri)) - y2 = yvert(vertex(2,itri)) - y3 = yvert(vertex(3,itri)) - lnod=0 - f3=-delta - do 30 i=1,iorder - f3 = f3+delta - f2 = -delta - f1 = 1.-f3+delta - do 20 j=i,iorder - lnod = lnod+1 - nod = node(lnod,itri) - f1 = f1-delta - f2 = f2+delta - if (xnode(nod).eq.-999.) then - xnode(nod) = f1*x1 + f2*x2 + f3*x3 - ynode(nod) = f1*y1 + f2*y2 + f3*y3 - endif - 20 continue - 30 continue - 40 continue -c -c sort to go first vertically then horizontally. use ipvt1 (not needed -c until factoring is done) to hold the inverse of the ordering vector -c simple bubble sort, @future use a better sort -c - eps = 4.*r1mach(4) - 1 continue - change = 0 - do 50 i=1,nnode0-1 - if (abs(xnode(ipvt1(i))-xnode(ipvt1(i+1))).lt.eps) then - if (ynode(ipvt1(i)).gt.ynode(ipvt1(i+1))) then - itemp = ipvt1(i) - ipvt1(i)=ipvt1(i+1) - ipvt1(i+1) = itemp - change=1 - endif - elseif (xnode(ipvt1(i)).gt.xnode(ipvt1(i+1))) then - itemp = ipvt1(i) - ipvt1(i)=ipvt1(i+1) - ipvt1(i+1) = itemp - change=1 - endif - 50 continue - if (change.eq.1) go to 1 -c -c set the ordering vector l1ord from its inverse in ipvt1 -c - do 60 i=1,nnode0 - l1ord(ipvt1(i))=i - 60 continue -c - if (outlev.ge.4) write(ioutpt,106) - 106 format(' determine bandwidth') -c -c determine bandwidth -c - do 70 i=1,nnode0 - do 70 j=1,mxidlo - id = idcoef(j,i) - if (id.ne.0) then - l1id=l1ord(id) - l1i =l1ord(i) - if (iabs(l1id-l1i) .gt. nband) nband=iabs(l1id-l1i) - endif - 70 continue -c -c verify enough space is allocated for coarse grid matrix bandwidth -c - if (nband .gt. ndband) then - write(ioutpt,100) ndband,nband - 100 format(/'FATAL ERROR -- not enough room for coarsest grid', - 1 ' factors for exact solve'/ - 2 'Current allocation is ndband = ',i6/ - 3 'Requires ndband >= ',i6/ - 4 'Change parameter statement in commons') - stop - endif -c - if (outlev.ge.3) write(ioutpt,108) - 108 format(' exact solve initializations complete') - return - end -c -c -------- exsolv -c - subroutine exsolv - include 'commons' -c -c solve exactly on the coarsest grid -c -c debug output -c - if (outlev.ge.3) write(ioutpt,100) - 100 format(' solve exactly on coarsest grid') -c -c exact solve on coarse grid by linpack band routine -c - if (outlev.ge.4) write(ioutpt,107) - 107 format(' copy coef to band form') -c -c initialize coef to 0. -c - do 80 j=1,nnode0 - do 80 i=1,3*nband+1 - coefl1(i,j)=0. - 80 continue -c -c copy coef to band form -c - do 95 i=1,nnode0 - l1i = l1ord(i) - if (bndnod(i)) then - coefl1(2*nband+1,l1i)=1. - else - coefl1(2*nband+1,l1i)=coef(mxidlo+1,i) - endif - do 90 j=1,mxidlo - id = idcoef(j,i) - if (id .ne. 0.) then - l1id = l1ord(id) - if (.not.bndnod(i)) - 1 coefl1(2*nband+1+l1i-l1id,l1id)=coef(j,i) - if (.not.bndnod(id)) - 1 coefl1(2*nband+1+l1id-l1i,l1i) =coef(j,i) - endif - 90 continue - 95 continue -c -c debug output matrix in band form -c - if (outlev.ge.5) then - write(ioutpt,184) - do 181 j=1,nnode0 - do 180 i=1,3*nband+1 - write(ioutpt,182) i,j,coefl1(i,j) - 180 continue - write(ioutpt,183) - 181 continue - 182 format(2i6,1pe15.8) - 183 format(' ') - 184 format(' matrix in band form:') - endif -c -c call linpack factorization -c - if (outlev.ge.4) write(ioutpt,102) - 102 format(' begin factoring coarse grid equations') -c - call sgbco(coefl1,3*ndband+1,nnode0,nband,nband,ipvt1,rcond,rs1) - if (rcond .lt. 10.*r1mach(4)) then - write(ioutpt,185) - endif - 185 format(/' WARNING -- The condition number of the coarse grid'/, - . ' matrix is very large.'/, - . ' The solution may be unreliable.'//) -c -c copy right side -c - do 10 i=1,nnode0 - if (bndnod(i)) then - rs1(l1ord(i)) = u(i) - else - rs1(l1ord(i)) = rs(i) - endif - 10 continue -c -c call linpack sgbsl -c - call sgbsl(coefl1,3*ndband+1,nnode0,nband,nband,ipvt1,rs1,0) -c -c copy solution -c - do 20 i=1,nnode0 - if (.not.bndnod(i)) u(i) = rs1(l1ord(i)) - 20 continue -c -c debug output -c - if (outlev.ge.5) then - write(ioutpt,200) - 200 format(' coarse grid solution is:') - call prvec(u,1,1) - endif -c - return - end -c -c -------- relax -c - subroutine relax(lev,nu) - include 'commons' - integer limit,lev,nu,vert,olnd,i,j,lim,blkhed -c -c regular form of relaxation -c -c perform nu half steps of red-black gauss siedel on level lev -c -c debug output -c - if (outlev.ge.4) write(ioutpt,901) nu,lev - 901 format(' relax ',i2,'/2 times on level ',i4) -c -c repeat nu/2 times -c - limit = int(nu/2) - if (limit.eq.0) go to 399 -c - do 320 i=1,limit -c - blkhed=-1 -c -c relax on highest level (red vertices) -c - vert = frstvt(lev) -c - 1 if (vert.eq.-1) go to 2 -c -c relax at the nodes associated with this vertex -c - call rlxred(vert) -c -c make note of which black nodes are neighbors -c - if (vrtlev(vert).gt.0) then - lim=iorder*iorder - else - lim=(iorder*(iorder+1))/2 - endif - do 20 j=1,lim - olnd = olndvt(j,vert) - if ((.not. inuse(olnd))) then - inuse(olnd) = .true. - nxtblk(olnd) = blkhed - blkhed=olnd - endif - 20 continue -c - vert = nextvt(vert) - go to 1 - 2 continue -c -c relax on other levels (black nodes) -c - olnd = blkhed - 101 if (olnd.eq.-1) go to 102 -c - call rlxblk(olnd) - inuse(olnd) = .false. - olnd = nxtblk(olnd) - go to 101 - 102 continue -c - 320 continue -c -c if nu is odd, relax red vertices one more time -c - 399 continue - if (2*limit .ne. nu) then - vert = frstvt(lev) - 401 if (vert.eq.-1) go to 402 -c -c relax at the nodes associated with this vertex -c - call rlxred(vert) -c - vert = nextvt(vert) - go to 401 - 402 continue - endif -c - return - end -c -c -------- rlxred -c - subroutine rlxred(vert) - include 'commons' - integer nod,i,limit,node1,k,vert - integer kpvt(ndord1),info - real sum -c -c relax at the red nodes associated with vert -c -c special case for iorder=2 -c - if (iorder.eq.2) then - if (.not.bndnod(vert)) then - sum=rs(vert) - do 5 i=1,4 - if (idcoef(i,vert).ne.0) then - sum=sum-coef(i,vert)*u(idcoef(i,vert)) - endif - 5 continue - u(vert)=sum/coef(5,vert) -c - endif - return - endif -c -c normal case for iorder>2 -c -c determine size of system -c - if (vrtlev(vert).gt.0) then - limit=nnodev - else - limit=nnodvb - endif - node1=nwndvt(vert) -c -c set up diagonal block associated with this vertex -c - do 110 i=1,limit - do 110 j=1,limit - block(i,j)=0. - 110 continue - do 130 i=1,limit - nod=node1+i-1 - do 120 j=1,mxidlo - if (idcoef(j,nod).ge.node1) then - block(i,idcoef(j,nod)-node1+1)=coef(j,nod) - block(idcoef(j,nod)-node1+1,i)=coef(j,nod) - endif - 120 continue - block(i,i)=coef(mxidlo+1,nod) - 130 continue -c -c set up right side for the system -c - do 20 k=1,limit -c - nod=node1-1+k -c - if (bndnod(nod)) go to 20 - u(nod)=rs(nod) -c -c contributions from lower triangular part -c - i=1 - 10 if (idcoef(i,nod).eq.0.or.idcoef(i,nod).ge.node1) go to 11 - u(nod)=u(nod)-coef(i,nod)*u(idcoef(i,nod)) - i=i+1 - go to 10 - 11 continue -c - 20 continue -c -c modify boundary rows -c - if (limit.eq.nnodvb) then - do 160 i=1,nnodvb - nod=i+node1-1 - if (bndnod(nod)) then - do 140 j=1,nnodvb - block(i,j)=0. - nod2=j+node1-1 - if (.not. bndnod(nod2)) then - u(nod2)=u(nod2)-block(j,i)*u(nod) - endif - block(j,i)=0. - 140 continue - block(i,i)=1. - endif - 160 continue - endif -c -c solve system -c - call ssifa(block,ndord1,limit,kpvt,info) - if (info .ne. 0) then - write(ioutpt,900) info - 900 format('ERROR -- linpack routine ssifa returned info = ',i8) - stop - endif - call ssisl(block,ndord1,limit,kpvt,u(node1)) -c - return - end -c -c -------- rlxblk -c - subroutine rlxblk(nod) - include 'commons' - integer nod,i,irow,icol - real sum -c -c relax at a black node -c - if (bndnod(nod)) return - sum=rs(nod) -c -c contributions from lower triangular part -c - do 10 i=1,mxidlo - if (idcoef(i,nod).ne.0) then - sum=sum-coef(i,nod)*u(idcoef(i,nod)) - endif - 10 continue -c -c contributions from upper triangular part -c - do 20 i=mxidlo+1,mxidup - if (idcoef(i,nod).ne.0) then -c -c find coef in lower triangular part -c - irow=idcoef(i,nod) - icol=1 - 1 if(idcoef(icol,irow).eq.nod) go to 2 - icol=icol+1 - if (icol.gt.mxidlo) then - write(ioutpt,100) nod,irow - stop - 100 format(/'********FATAL ERROR********'// - 1 ' couldnt find column ',i10,' in row ',i10, - 2 ' in rlxblk') - endif - go to 1 - 2 continue - sum=sum-coef(icol,irow)*u(irow) - endif - 20 continue -c - u(nod)=sum/coef(mxidlo+1,nod) -c - return - end -c -c -------- locrlx -c - subroutine locrlx(t) - include 'commons' - integer i,j,n,tri,vrt,limit,nod,t,trlist(24),neftri - logical onlist -c -c local relaxations after adding a new vertex -c the relaxation at the new nodes has already been done -c relaxes at the old nodes in the new triangles -c -c debug output -c - if (outlev.ge.4) write(ioutpt,100) - 100 format(' local relaxations') -c -c limit=number of old nodes -c - if (vrtlev(nvert).lt.0) then - limit=(iorder*(iorder+1))/2 - else - limit=iorder*iorder - endif -c -c relax at old nodes -c - do 10 n=1,limit - nod=olndvt(n,nvert) - call rlxblk(nod) - 10 continue -c -c correct error indicators at triangles by old nodes -c - if (vrtlev(nvert).lt.0) then - limit=3 - else - limit=4 - endif -c -c make a list of effected triangles -c - neftri=0 - do 30 i=1,limit -c -c find vertex -c - if (i.eq.1) then - vrt=vertex(1,t) - elseif (i.eq.2) then - vrt=vertex(2,t) - elseif (i.eq.3) then - vrt=vertex(1,ntri) - else - vrt=vertex(2,ntri) - endif -c -c add neighboring triangles to list if not already there -c - do 20 j=1,8 - tri=tringl(j,vrt) - if (tri.eq.0) go to 20 - onlist=.false. - if (neftri.ne.0) then - do 15 k=1,neftri - if (trlist(k).eq.tri) onlist=.true. - 15 continue - endif - if (.not. onlist) then - neftri=neftri+1 - trlist(neftri)=tri - endif - 20 continue - 30 continue -c -c debug output -c - if (outlev.ge.5) write(ioutpt,200) (trlist(k),k=1,neftri) - 200 format(' triangles whose error indicator is effected ', - 1 'by local relaxation '/12i6/12i6) -c -c correct error indicator for each neighboring triangle -c - do 40 j=1,neftri - tri=trlist(j) - call eitri(tri) - if(coefei(iorder*iorder+1,1,tri).eq.0.) - 1 tri=neigh(3,tri) - call elstrm(tri) - call elstad(tri) - 40 continue - return - end -c -c -------- esterr -c - subroutine esterr - include 'commons' - integer t,i - external clock -c -c compute error indicator for each triangle using -c 1 point local dirichlet, and compute global error estimate -c -c debug output -c - if (outlev.ge.2) write(ioutpt,101) - 101 format(/' begin error indicators') - if (outlev.ge.4) write(ioutpt,102) - 102 format(/' error indicators are:') -c -c begin timing -c - timeel = clock(0) -c -c compute error indicators -c - eimax = 0. - do 10 t=1,ntri - if (coefei(iorder*iorder+1,1,t).ne.0.) then - call eitri(t) - if (errind(t).gt.eimax) eimax=errind(t) - endif -c -c debug output -c - if (outlev.ge.4) write(ioutpt,103) t,errind(t) - 103 format(1x,i6,2x,1pe15.8) -c - 10 continue -c -c compute global error estimate -c - gerest = 0. - do 18 t=1,ntri - if (coefei(iorder*iorder+1,1,t).ne.0.) then - if (neigh(3,t).le.0) then - gerest=gerest+errind(t)*errind(t) - elseif(neigh(3,neigh(3,t)).eq.t) then - gerest=gerest+errind(t)*errind(t) - else - gerest=gerest+errind(t)*errind(t)/2. - endif - endif - 18 continue -c -c global error estimates -c - gerest = sqrt(abs((2**(iorder-1))*gerest/(2.**(iorder-1)-1.))) - if (abs(unrm) .gt. 1.e-10) then - rerest = gerest/unrm - else - rerest = gerest - endif -c -c create lists of error indicators in various ranges -c -c empty lists -c - do 20 i=1,4 - eihead(i) = -1 - eitail(i) = -1 - 20 continue -c -c add triangles with error indicators to lists -c - do 30 t=1,ntri - if (coefei(iorder*iorder+1,1,t).eq.0.) then - nxttri(t) = -1 - pretri(t) = -1 - else - call elstad(t) - endif - 30 continue -c -c finish timing -c - timeel = clock(0) - timeel - timeet = timeet + timeel -c -c save error estimate for gnuplot convergence file -c - gpeest(gplev)=rerest - gptime(gplev)=clock(0)-timett - gplev=gplev+1 -c -c debug output -c - if (outlev.ge.4) then - write(ioutpt,111) - do 110 i=1,4 - write(ioutpt,112) i - t = eihead(i) - 121 if (t.eq.-1) go to 122 - write(ioutpt,113) t - t = nxttri(t) - go to 121 - 122 continue - 110 continue - endif - 111 format(/' error indicator lists are') - 112 format(' list ',i2) - 113 format(i6) -c - if (outlev.ge.2) then - write(ioutpt,104) - write(ioutpt,105) eimax,gerest - if (abs(gerr) .gt. 1.e-10) then - effind = gerest/gerr - else - effind = 1. - endif - if (abs(rgerr) .gt. 1.e-10) then - refind = rerest/rgerr - else - refind = 1. - endif - write(ioutpt,106) effind - write(ioutpt,116) rerest,refind - write(ioutpt,107) timeel,timeet - endif - if (outlev.ge.4) then - call plttri(0,3) - endif - 104 format(/' error indicators and estimates complete') - 105 format(/' maximum error indicator ',1pe15.8/ - 1 ' error estimate ',1pe15.8) - 106 format( ' effectivity index ',1pe15.8) - 116 format( ' relative error estimate ',1pe15.8/ - 1 ' relative effect index ',1pe15.8) - 107 format(/' time for error estimates (this grid) ',f10.2/ - 1 ' time for error estimates (all grids) ',f10.2) -c -c write data files for run time graphics, check the menu, and -c update the plots -c - if (outlev.ge.2 .and. grafic) write(ioutpt,200) - 200 format(' begin update to convergence graphics') - if (grafic) then - call filcon - if (menuon) call mn2mg - call grfcon - endif - if (outlev.ge.2 .and. grafic) write(ioutpt,300) - 300 format(' graphics complete') -c - return - end -c -c -------- eitri -c - subroutine eitri(t) - include 'commons' - integer t,t1,nbr,ord2,nbrtyp,i,nod,j,limit,id,sumord,bctype - integer kpvt(ndord1) - real sum,bcrhs,cu - logical bondry,cmpdiv,dirich -c -c compute error indicator for triangle t -c -c avoid null triangle - if (t.eq.0) return -c -c determine whether t or its mate has the problem -c - ord2 = iorder*iorder - if (coefei(ord2+1,1,t).eq.0.) then - t1=neigh(3,t) - else - t1=t - endif -c -c set neighbor and case flags -c - nbr=neigh(3,t1) - if (nbr.le.0) then - bondry=.true. -c@@@ call is just to get bctype - call bcond(xvert(vertex(1,t1)),yvert(vertex(1,t1)), - 1 -nbr,cu,bcrhs,bctype) - if (bctype.eq.1) then - dirich = .true. - else - dirich = .false. - endif - cmpdiv=.true. - nbrtyp=0 - else - dirich = .false. - bondry=.false. - if (neigh(3,nbr).eq.t1) then - cmpdiv=.true. - nbrtyp=0 - else - cmpdiv=.false. - if (vertex(1,nbr).eq.vertex(1,t1)) then - nbrtyp=1 - else - nbrtyp=2 - endif - endif - endif -c -c special case for iorder=2 -c - if (iorder.eq.2) then -c -c boundary triangle - if (dirich) then - uei(1)=bcei(1,t1)-(u(vertex(1,t1))+u(vertex(2,t1)))/2. -c interior triangle - else -c contribution from nodes in t1 - uei(1)=rsei(1,t1)-coefei(1,1,t1)*u(vertex(1,t1)) - 1 -coefei(2,1,t1)*u(vertex(2,t1)) - 2 -coefei(3,1,t1)*u(vertex(3,t1)) -c contribution from node in nbr depends on compatability -c (and aren't there if the edge is boundary) - if (.not.bondry) then - if (cmpdiv) then - uei(1)=uei(1)-coefei(4,1,t1)*u(vertex(3,nbr)) - else - uei(1)=uei(1)-coefei(4,1,t1)*(u(vertex(1,nbr))+ - 1 u(vertex(2,nbr)))/2. - endif - endif -c solve system - uei(1)=uei(1)/coefei(5,1,t1) -c convert to coefficent of 2-level basis - uei(1)=uei(1)-(u(vertex(1,t1))+u(vertex(2,t1)))/2. - endif -c -c compute energy norm -c - errind(t1)=sqrt(abs(uei(1)*uei(1)*coefei(5,1,t1))) -c - else -c -c higher order methods -c -c set solution values for the neighbor if not compatably divisible -c - if (.not. cmpdiv) then - do 20 i=1,(iorder*(iorder-1))/2 -c -c find local numbering of node -c - nod=renum(nbrtyp,i+iorder) - if (nod.gt.0) then -c node is an existing node in the neighbor - ueiold(i)=u(node(nod,nbr)) - else -c node will be added when the neighbor is divided -c evaluate solution there - nod=-nod - ueiold(i)=0. - do 10 j=1,nbasch(nod) - ueiold(i)=ueiold(i)+cbasch(j,nod) - 1 *u(node(ibasch(j,nod),nbr)) - 10 continue - endif - 20 continue - endif -c -c set up right side for system -c - if (bondry) then - limit=nnodvb - else - limit=nnodev - endif - do 40 i=1,limit - uei(i)=rsei(i,t1)+rsei2(i,t1) -c -c no changes for nodes on the boundary -c - if (.not. bondry .or. nbasch(i).ne.iorder - 1 .or. .not. dirich) then -c -c contributions from existing nodes -c - j=1 - 21 if (j.gt.ord2) go to 22 - if (idcoei(j,i,t1).eq.0) go to 22 - id=idcoei(j,i,t1) - if (id.lt.0) then -c node is in non compatable neighbor - uei(i)=uei(i)-coefei(j,i,t1)*ueiold(-id) - else -c node exists - uei(i)=uei(i)-coefei(j,i,t1)*u(id) - endif - j=j+1 - go to 21 - 22 continue - else -c -c set Dirichlet boundary condition -c - uei(i)=bcei(i,t1) - endif - 40 continue -c -c copy blokei to block -c - do 51 i=1,limit - do 51 j=1,limit - block(i,j)=blokei(i,j,t1) - 51 continue -c -c solve system -c - call ssifa(block,ndord1,limit,kpvt,info) - if (info .ne. 0) then - write(ioutpt,900) info - 900 format('ERROR -- linpack routine ssifa returned info = ',i8) - stop - endif - call ssisl(block,ndord1,limit,kpvt,uei) -c -c convert solution to 2-level basis coefficents -c - sumord=(iorder*(iorder+1))/2 - do 120 i=1,limit - do 110 j=1,nbasch(i) - if (ibasch(j,i).le.sumord) then - uei(i)=uei(i)-cbasch(j,i)*u(node(ibasch(j,i),t1)) - elseif (cmpdiv) then - uei(i)=uei(i)-cbasch(j,i) - 1 *u(node(ibasch(j,i)-sumord+iorder,nbr)) - else - uei(i)=uei(i)-cbasch(j,i)*ueiold(ibasch(j,i)-sumord) - endif - 110 continue - 120 continue -c -c compute energy norm -c - errind(t1)=0. - sum=0. - do 140 i=1,limit - do 130 j=1,limit - sum=sum+blokei(i,j,t1)*uei(i)*uei(j) - 130 continue - 140 continue - errind(t1)=sqrt(abs(sum)) -c - endif -c - return - end -c -c -------- elstad -c - subroutine elstad(t) - include 'commons' - integer t,list - real e,cut1,cut2,cut3 - logical top -c - if (unifrm) return -c -c add triangle t to the correct error indicator list -c - e = errind(t) -c -c if e>eimax, shift lists to create a new list of largest error -c indicators. Don't worry about needing to shift more than once -c @future examine the effect of only checking for shift once -c - if ( e.gt.eimax ) then -c combine last two lists - if (eihead(3).ne.-1) then - if (eihead(4).eq.-1) then - eihead(4)=eihead(3) - eitail(4)=eitail(3) - else - nxttri(eitail(3)) = eihead(4) - pretri(eihead(4)) = eitail(3) - eihead(4) = eihead(3) - endif - endif -c shift other lists - eihead(3) = eihead(2) - eitail(3) = eitail(2) - eihead(2) = eihead(1) - eitail(2) = eitail(1) - eihead(1) = -1 - eitail(1) = -1 -c increase eimax - eimax = eimax*mgfreq**(iorder/2.) - if (outlev.ge.4) write(ioutpt,100) - 100 format(' shift error indicator lists for larger eimax') - endif -c -c determine which list this triangle goes in and whether -c it goes in the top or the bottom -c - cut1 = eimax*mgfreq**(-iorder/2.) - cut2 = eimax*mgfreq**(-float(iorder)) - cut3 = eimax*mgfreq**(-3.*iorder/2.) - if (e.ge.cut1) then - list = 1 - top = (e.gt. (eimax+cut1)/2.) - elseif (e.ge.cut2) then - list=2 - top = (e.gt. (cut1+cut2)/2.) - elseif (e.ge.cut3) then - list=3 - top = (e.gt. (cut2+cut3)/2.) - else - list=4 - top = (e.gt. cut3/2.) - endif -c -c add triangle to list -c - if (top) then -c - nxttri(t) = eihead(list) - pretri(t) = -1 - if (eihead(list) .eq. -1) then - eitail(list) = t - else - pretri(eihead(list)) = t - endif - eihead(list) = t -c - else -c - pretri(t) = eitail(list) - nxttri(t) = -1 - if (eitail(list) .eq. -1) then - eihead(list) = t - else - nxttri(eitail(list)) = t - endif - eitail(list) = t -c - endif -c -c debug output -c - if (outlev .ge. 5) then - if (top) then - write(ioutpt,201) t,list - else - write(ioutpt,202) t,list - endif - endif - 201 format(' add triangle ',i6,' to top of error indicator', - 1 ' list ',i2) - 202 format(' add triangle ',i6,' to bottom of error indicator', - 1 ' list ',i2) -c - return - end -c -c -------- elstrm -c - subroutine elstrm(t) - include 'commons' - integer t,list -c - if (unifrm) return -c -c remove triangle t from the error indicator lists -c -c if t is at the head of a list, find out which list -c - if (pretri(t).eq.-1) then - if (eihead(1).eq.t) then - list=1 - elseif (eihead(2).eq.t) then - list=2 - elseif (eihead(3).eq.t) then - list=3 - elseif (eihead(4).eq.t) then - list=4 - else - write(ioutpt,100) t - stop - 100 format(' ********FATAL ERROR********'// - 1 ' couldnt find which list ',i6,' is the head of') - endif -c -c change pointer before t -c - eihead(list) = nxttri(t) - else - nxttri(pretri(t)) = nxttri(t) - endif -c -c if t is at the end of a list, find out which list -c - if (nxttri(t) .eq. -1) then - if (eitail(1) .eq. t) then - list = 1 - elseif (eitail(2) .eq. t) then - list = 2 - elseif (eitail(3) .eq. t) then - list = 3 - elseif (eitail(4) .eq. t) then - list = 4 - else - write(ioutpt,200) t - stop - 200 format(' ********FATAL ERROR********'// - 1 ' couldnt find which list ',i6,' is the tail of') - endif -c -c change back pointer after t -c - eitail(list) = pretri(t) - else - pretri(nxttri(t)) = pretri(t) - endif -c -c debug output -c - if (outlev.ge.5) write(ioutpt,300) t - 300 format(' remove triangle ',i6,' from error indicator lists') -c - return - end -c -c -------- errors -c - subroutine errors(emaxv,emaxn,emaxq,eenrg) - include 'commons' - integer t,denom,i1,i2,i3,i,nod,j - real emaxv,emaxn,emaxq,eenrg,x1,x2,x3,y1,y2,y3,area, - 1 det,dz1dx,dz2dx,dz3dx,dz1dy,dz2dy,dz3dy,x,y - real uval, uxval, uyval, err, errx, erry - real qpf,qpp,qpq,qpr -c -c compute norms of error -c 4 norms : max norm at vertices, nodes, and quadrature points -c and continuous energy norm computed with the same -c quadrature rule as used for computing inner products -c - emaxv=0. - emaxn=0. - emaxq=0. - eenrg=0. -c -c computations one triangle at a time -c - do 40 t=1,ntri -c -c x-y coordinates of vertices of this triangle -c - x1=xvert(vertex(1,t)) - x2=xvert(vertex(2,t)) - x3=xvert(vertex(3,t)) - y1=yvert(vertex(1,t)) - y2=yvert(vertex(2,t)) - y3=yvert(vertex(3,t)) -c -c area of triangle and derivatives of zeta functions wrt x and y -c - det = x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2) - area = abs(det/2.) - dz1dx = (y2-y3)/det - dz2dx = (y3-y1)/det - dz3dx = (y1-y2)/det - dz1dy = (x3-x2)/det - dz2dy = (x1-x3)/det - dz3dy = (x2-x1)/det -c -c weights for computing x-y coordinates of nodes -c - denom = iorder-1 - i1 = denom+1 - i2 = -1 - i3 = 0 -c -c pass through nodes for max err at nodes and vertices -c - do 10 i=1,nnodtr - nod = node(i,t) -c -c x-y coordinates of node -c - i1=i1-1 - i2=i2+1 - if (i1.lt.0) then - i3=i3+1 - i1=denom-i3 - i2=0 - endif - x=(i1*x1 + i2*x2 + i3*x3)/float(denom) - y=(i1*y1 + i2*y2 + i3*y3)/float(denom) -c -c error at this node -c - err = u(nod)-true(x,y) -c -c check for max -c - if (abs(err).gt.emaxn) emaxn=abs(err) - if ((i.eq.1 .or. i.eq.iorder .or. i.eq.nnodtr) - 1 .and. abs(err).gt.emaxv) emaxv=abs(err) -c - 10 continue -c -c pass through quadrature points for max err at quad pts and -c to compute integral for energy norm -c - do 30 i=1,nqpt -c -c x-y coordinates of quadrature points -c - x=x1*quadpt(1,i)+x2*quadpt(2,i)+x3*quadpt(3,i) - y=y1*quadpt(1,i)+y2*quadpt(2,i)+y3*quadpt(3,i) -c -c pde coefficents at quadrature point -c - call pde(x,y,qpp,qpq,qpr,qpf) -c -c value of approximate solution and derivatives at quadrature point -c - uval = 0. - uxval= 0. - uyval= 0. - do 20 j=1,nnodtr - uval = uval+qpbas(j,i)*u(node(j,t)) - uxval= uxval+(qpdbdz(1,j,i)*dz1dx - 1 + qpdbdz(2,j,i)*dz2dx - 2 + qpdbdz(3,j,i)*dz3dx)*u(node(j,t)) - uyval= uyval+(qpdbdz(1,j,i)*dz1dy - 1 + qpdbdz(2,j,i)*dz2dy - 2 + qpdbdz(3,j,i)*dz3dy)*u(node(j,t)) - 20 continue -c -c error and derivatives of error at quadrature point -c - err = uval - true(x,y) - errx = uxval - truex(x,y) - erry = uyval - truey(x,y) -c -c check for max error -c - if (abs(err).gt.emaxq) emaxq=abs(err) -c -c contribution to integral -c - eenrg = eenrg + area*quadw(i)* - 1 abs(qpp*errx*errx + qpq*erry*erry - 2 + qpr*err*err) -c - 30 continue - 40 continue -c - eenrg = sqrt(eenrg) -c - return - end -c -c -------- errred -c - subroutine errred(loop,oemaxv,oemaxn,oemaxq,oeenrg) - include 'commons' - integer loop - real emaxv,oemaxv,remaxv,emaxn,oemaxn,remax