********************************************************************************* * * * Grid optimization code * * Developed by Joseph Yinglong Zhang * * Supervised by Prof. A.M. Baptista * * Center for Coastal and Land-Margin Research * * Oregon Graduate Institute of Science and Technology * * Version 5b, Nov 2002 * * * * AGrid is a 2D grid optimization code that seeks to * * equi-distribute edge errors, given an error field associated * * with a grid. It utilizes five core modules to achieve the * * error equi-distribution: move node, refine, coarse and nswap edge, * * and snap node to an edge. * * * * * ********************************************************************************* ********************************************************************************* * * * Version updates (from v.5): * * (1) Split the 'coarsen' part in coarsen_geo into 2: eliminating * * 3-node balls and coarsen; added virus scan after each part * * for easier debugging. * * (2) Change the arguments of tridag and cyclic for robustness. * * (3) Changed input/output file names to more meaningful ones. * * * ********************************************************************************* * Input files: * (1) bg grid file (bg.gr3): standard grid file in gredit format, * with the depth being the error at the node; * (2) fg grid file (fg.gr3): standard grid file in gredit format, * depth not used; * (3) boundary nodes file (fg.ebn): in gredit format (edit boundary nodes); * (4) main parameter input file (param.in). * (5) reserved nodes list (rsv_nd.in with irsvflag=2): * (6) reserved elements list (rsv_elem.in ): element property format (like dpe.dat * in Elcirc). * Output files: * (1) screen output & main output file (fort.10); * (2) edge error distribution file (err_dis.dat); * (3) intermediate grid and output files (tmp.gr3 and store.dat); * (4) optimized grid (op.gr3); * (5) optimized grid after geometric correction (final.gr3 and final.ll); * (6) reserved nodes (rsv_nd.out); * (7) reserved element list for final grid (rsv_elem.out), if rsv_elem.in * is part of input (i.e. ires_elem=1). * param.in (paramter file): * ALNUM= warning message; * ics,slam0,sfea0: coordinate system. ics=1: cartesian; =2: lat long, * and slam0, sfea0 are center of projection * (in degrees). For the latter option, there are 2 * output files, converged.gr being the cartesian * file, and converged.ll being the lat long file. * ihot: flag for hot start; =0: no hot start; =1: hot start (can only * start from the first optimiation part, i.e., not in geometry part). * icoflag emin emax: flag; =0: error cutoff disabled; =1: enabled, * and emin and emax are lower and upper error * cutoff values; * ismflag it_sm: flag; =0: no error field smoothing; =1: error * smoothing using least square method and it_sm is the * no. of iterations in smoothing (ignored when * ismflag=0); * ibg: =0: fort.14 identical to fort.13; =0: read in fort.14. * cst_neg: smallest legal area (in lieu of 0). * irsvflag: =0: no reserved nodes; =1: reserve based on curvature; * =2: reserve according to fort.17. * if(irsvflag=1), degrm: maximum change of slope (in degrees) * allowed in each boundary interval. * ires_elem: =0: no reserved elements; =1: reserve elements according to fort.19. * ibifur: flag for automatic search for thresholds rlam1-3. =0:disabled; * =1: enabled. * if ibifur=0, next line will be: * icsflag rlam1-3: flag; =0: relative thresholds for refine, * coarsen and snap * modules; =1: absolute thresholds; =2: absolute * thresholds as ratios to INITIAL averaged error. * rlam1-3: thresholds in refine, coarsen and snap modules; * if ibifur=1, next line will be: * ndmin ndmax: target interval for the # of nodes. * iconflag: flag for convergence criterion used in optimization; * =0: absolute convergence; =1: relaxed convergence. * if(iconflag=1), dn_r,an_r,se_r: the three tolerances * give the maximum ratios of deleted nodes, * added nodes and swapped edges allowed for * convergence. * iflaggeo: flag for geometric optimization to get rid of very * skew elements; =0: disabled (then the remaining * parameters are all ignored); =1: enabled. * askew: acceptable skewness (skewness=longest edge/sqrt(area/pi)). * it_geo: no. of iterations for geometric run. * Global arrays: xn, yn, ic1, irsve, nrsv, ibdf, itab, nsur_e, ic_e *********************************************************************************** * program AGrid include 'agrid.h' dimension nbd(mnp),ndiv(mnb),error1(mnp),iwork(10) C... Following matrices used in calculating curvature dimension aa(mnp),bb(mnp),cc(mnp),rrx(mnp),rry(mnp) dimension ddx(mnp),ddy(mnp),curv(mnp) C... Initialize data ibdf/mnp*0/ data nrsv/mnp*0/ data irsve/mne*0/ data iskew/mne*0/ C... Open parameter file (fort.15) C... open(15,file='param.in',status='old') read(15,*) read(15,*)ics,slam0,sfea0 if(ics.ne.1.and.ics.ne.2) then write(*,*)'Unknown ics',ics stop endif slam0=slam0*pi/180 sfea0=sfea0*pi/180 rearth=6378206.4d0 read(15,*)ihot if(ihot.ne.0.and.ihot.ne.1) then write(*,*)'Unknown ihot',ihot stop endif read(15,*)icoflag,emin,emax if(icoflag.ne.0.and.icoflag.ne.1) then write(*,*)'Unknown icoflag',icoflag stop endif read(15,*)ismflag,it_sm if(ismflag.ne.0.and.ismflag.ne.1) then write(*,*)'Unknown ismflag',ismflag stop endif read(15,*)ibg if(ibg.ne.0.and.ibg.ne.1) then write(*,*)'Unknown ibg',ibg stop endif read(15,*)cst_neg read(15,*)irsvflag if(irsvflag.lt.0.or.irsvflag.gt.2) then write(*,*)'Unknown irsvflag',irsvflag stop endif if(irsvflag.eq.1) read(15,*)degrm read(15,*)ires_elem if(ires_elem.ne.0.and.ires_elem.ne.1) then write(*,*)'Unknown ires_elem',ires_elem stop endif read(15,*)ibifur if(ibifur.eq.0) then read(15,*)icsflag,rlam1,rlam2,rlam3 if(icsflag.lt.0.or.icsflag.gt.2) then write(*,*)'Unknown icsflag',icsflag stop endif else if(ibifur.eq.1) then read(15,*)ndmin,ndmax,rat_u,rat_l if(ndmax.gt.mnp) then write(*,*)'maximum # of nodes exceeded' stop endif if(ndmin.gt.ndmax) then write(*,*)'ndmin > ndmax' stop endif if(rat_u.le.1.or.rat_l.le.0.or.rat_l.ge.1) then write(*,*)'rat_u must be > 1, 0 < rat_l < 1' stop endif else write(*,*)'Unknown ibifur' stop endif read(15,*)iconflag if(iconflag.ne.0.and.iconflag.ne.1) then write(*,*)'Unknown iconflag',iconflag stop endif if(iconflag.eq.1) read(15,*)dn_r,an_r,se_r read(15,*)iflaggeo if(iflaggeo.ne.0.and.iflaggeo.ne.1) then write(*,*)'Unknown iflaggeo',iflaggeo stop endif read(15,*)askew read(15,*)it_geo close(15) C... Read in bg grid with error field C... open(13,file='bg.gr3',status='old') read(13,*) read(13,*)nt0,nn0 if(nn0.gt.mnp.or.nt0.gt.mne) then write(*,*)'nn0 > mnp or nt0 > mne' stop endif do i=1,nn0 read(13,*)j,x0(i),y0(i),error0(i) C Convert to Cartesian coordinates if(ics.eq.2) then x0(i)=x0(i)*pi/180 y0(i)=y0(i)*pi/180 x0(i)=rearth*(x0(i)-slam0)*dcos(sfea0) y0(i)=rearth*y0(i) endif if(error0(i).lt.0.0d0) then write(*,*)'Negative error at node',i write(*,*)'Check fort.13!' stop endif if(icoflag.eq.1) error0(i)=dmin1(dmax1(error0(i),emin),emax) C Max. & min. if(i.eq.1) then xmin=x0(1) ymin=y0(1) xmax=x0(1) ymax=y0(1) else xmin=dmin1(xmin,x0(i)) xmax=dmax1(xmax,x0(i)) ymin=dmin1(ymin,y0(i)) ymax=dmax1(ymax,y0(i)) endif enddo !i=1,nn0 do i=1,nt0 read(13,*)j,k,(ic0(i,l),l=1,3) C negative elements n1=ic0(i,1) n2=ic0(i,2) n3=ic0(i,3) area0(i)=area(x0(n1),x0(n2),x0(n3),y0(n1),y0(n2),y0(n3)) if(area0(i).le.0) then write(*,*)'Negative element',i stop endif enddo close(13) C Error smoothing if(ismflag.ne.0) then call surrounding2 do it=1,it_sm call efilter(error1) do i=1,nn0 if(icoflag.eq.0) then error0(i)=error1(i) else error0(i)=dmin1(dmax1(error0(i),emin),emax) endif enddo !i enddo !it C Output initial errors open(12,file='init_bg.gr3') write(12,*)'Initial background grid with smoothed errors' write(12,*)nt0,nn0 do i=1,nn0 write(12,100)i,x0(i),y0(i),error0(i) enddo do i=1,nt0 write(12,101)i,3,(ic0(i,j),j=1,3) enddo close(12) endif !ismflag.ne.0 C Compute ic4 call surrounding(0) do i=1,nt0 do j=1,3 n1=ic0(i,j) n2=ic0(i,mod(j,3)+1) ic4(i,j)=0 !bnd flag icount=0 do l=1,nsur_e(n1) ne=ic_e(n1,l) if(ne.ne.i.and. & (ic0(ne,1).eq.n2.or.ic0(ne,2).eq.n2.or.ic0(ne,3).eq.n2)) then icount=icount+1 if(icount.gt.1) then write(*,*)'More than 2 elements sharing same edge in bg' write(*,*)n1,n2 stop endif ic4(i,j)=ne endif enddo !l=1,nsur_e(n1) enddo !j=1,3 enddo !i=1,nt0 C Cold and hot start if(ihot.eq.1) then C---------------------------------------------------------------------------- if(ibifur.ne.0) then write(*,*)'no automatic search for hot start' stop endif if(icsflag.eq.2) then write(*,*)'Please change the flag to 1 and use absolute' write(*,*)'values for rlam1-3 (see the beginning of' write(*,*)'fort.10 for their last-used values)' stop endif call readin(istep) C---------------------------------------------------------------------------- else !cold start istep=1 C... Read in edit grid C... if(ibg.eq.0) then nn=nn0 do i=1,nn xn(i)=x0(i) yn(i)=y0(i) enddo nt=nt0 do i=1,nt do j=1,3 ic1(i,j)=ic0(i,j) enddo enddo else !ibg/=0; read in fort.14 open(14,file='fg.gr3',status='old') read(14,*) read(14,*)nt,nn if(nn.gt.mnp.or.nt.gt.mne) then write(*,*)'nn > mnp or nt > mne' stop endif do i=1,nn read(14,*)j,xn(i),yn(i),dum C Convert to Cartrsian if(ics.eq.2) then xn(i)=xn(i)*pi/180 yn(i)=yn(i)*pi/180 xn(i)=rearth*(xn(i)-slam0)*dcos(sfea0) yn(i)=rearth*yn(i) endif enddo !i=1,nn do i=1,nt read(14,*)j,k,(ic1(i,l),l=1,3) C Negative elements n1=ic1(i,1) n2=ic1(i,2) n3=ic1(i,3) areai=area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)) if(areai.le.0) then write(*,*)'Negative element in edit grid',i stop endif enddo !i close(14) endif !ibg C Ball call surrounding(1) C Initialize search table if(ibg.eq.0) then do i=1,nn itab(i)=ic_e(i,1) enddo else C more work write(*,*)'Under construction' stop endif C... Read in boundary nodes open(16,file='fg.ebn',status='old') read(16,*) read(16,*)noofbd if(noofbd.gt.mnb) then write(*,*)'Maximum # of bnd exceeded!' stop endif icount=0 do i=1,noofbd read(16,*)ndiv(i) do j=1,ndiv(i) icount=icount+1 read(16,*)nbd(icount) ibdf(nbd(icount))=1 !flag enddo !j enddo !i close(16) nbd_m=icount C... Reserving bnd nodes C... if(irsvflag.eq.1) then C Reserve nodes based on curvature C Calculate curvature for the bnd nstart=0 do ibd=1,noofbd C Form matrix do i=1,ndiv(ibd) aa(i)=1 bb(i)=4 cc(i)=1 if(i.eq.1) then in1=nbd(ndiv(ibd)+nstart) in2=nbd(i+nstart) in3=nbd(i+1+nstart) else if(i.eq.ndiv(ibd)) then in1=nbd(i-1+nstart) in2=nbd(i+nstart) in3=nbd(1+nstart) else in1=nbd(i-1+nstart) in2=nbd(i+nstart) in3=nbd(i+1+nstart) endif rrx(i)=6*(xn(in3)+xn(in1)-2*xn(in2)) rry(i)=6*(yn(in3)+yn(in1)-2*yn(in2)) enddo !i=1,ndiv(ibd) alpha=1 beta=1 call cyclic(aa,bb,cc,alpha,beta,rrx,ddx,ndiv(ibd)) call cyclic(aa,bb,cc,alpha,beta,rry,ddy,ndiv(ibd)) curvmax=0 do i=1,ndiv(ibd) in2=nbd(i+nstart) if(i.eq.ndiv(ibd)) then in3=nbd(1+nstart) xdot=xn(in3)-xn(in2)-ddx(i)/3-ddx(1)/6 ydot=yn(in3)-yn(in2)-ddy(i)/3-ddy(1)/6 else in3=nbd(i+1+nstart) xdot=xn(in3)-xn(in2)-ddx(i)/3-ddx(i+1)/6 ydot=yn(in3)-yn(in2)-ddy(i)/3-ddy(i+1)/6 endif if(xdot*xdot+ydot*ydot.eq.0) then curv(i)=1.0d4 write(*,*)'Warning: large curvature at node',in2 else curv(i)=dabs(xdot*ddy(i)-ydot*ddx(i))/ & dsqrt(xdot*xdot+ydot*ydot)**3 endif if(curv(i).gt.curvmax) then curvmax=curv(i) i0=i endif enddo !i=1,ndiv(ibd) nrsv(nbd(i0+nstart))=1 i0i=i0-1 !left end pt cdeg=0 23 continue i0i=i0i+1 i0i1=i0i+1 !right end pt if(i0i.gt.ndiv(ibd)) i0i=i0i-ndiv(ibd) if(i0i1.gt.ndiv(ibd)) i0i1=i0i1-ndiv(ibd) C Check if i0i1 goes back to i0 if(i0i1.eq.i0) go to 24 in2=nbd(i0i+nstart) in3=nbd(i0i1+nstart) rl=dsqrt((xn(in3)-xn(in2))**2+(yn(in3)-yn(in2))**2) cdeg=cdeg+rl*0.5d0*(curv(i0i)+curv(i0i1)) if(cdeg.gt.degrm*pi/180) then nrsv(in3)=1 i0i=i0i1-1 cdeg=0 go to 23 endif go to 23 24 continue nstart=nstart+ndiv(ibd) enddo !ibd=1,noofbd C Output reserved nodes to rsv.dat (in build pts) open(17,file='rsv_nd.out') write(17,*)'Reserved nodes' write(17,*)nn do i=1,nn write(17,100)i,xn(i),yn(i),float(nrsv(i)) enddo close(17) endif !irsvflag.eq.1 C... Other ways of reserving nodes C... if(irsvflag.eq.0) write(10,*)'No reserved nodes!' if(irsvflag.eq.2) then !read in from rsv_nd.in open(17,file='rsv_nd.in',status='old') read(17,*)nrsv_m do i=1,nrsv_m read(17,*)nd nrsv(nd)=1 enddo close(17) endif C... Read in reserved element list C... if(ires_elem.eq.0) then write(*,*)'No elements reserved' else !ires_elem=1 open(19,file='rsv_elem.in',status='old') nrsve=0 do i=1,nt read(19,*)j,dum if(dum.gt.0.5) then irsve(i)=1 nrsve=nrsve+1 endif enddo write(*,*)nrsve,' elements reserved' close(19) endif C Virus scan C Check bnd do i=1,nbd_m do j=i+1,nbd_m if(nbd(i).eq.nbd(j)) then write(*,*)'Redundant bd nodes!',nbd(i) stop endif enddo !j enddo !i call vscan C... Check integrity of the bnd call integrity write(*,*)'Virus scan passed' C Compute initial edges errors distribution call err_distri2(eaver,sd,err_min,err_max) write(*,*)'Average edge error=',eaver write(10,*)'Average edge error=',eaver write(10,*)'Error bounds=',err_max,err_min if(err_min.ne.0) write(10,*)'Error ratio=',err_max/err_min C Auto search if(ibifur.eq.1) then write(*,*)'Automatic search for rlam1-3' write(10,*)'Automatic search for rlam1-3' icsflag=1 rlam1=2*eaver rlam2=0.5*eaver rlam3=0.4*eaver call store(istep) endif C Redefine rlam1-3 if(icsflag.eq.2) then rlam1=rlam1*eaver rlam2=rlam2*eaver rlam3=rlam3*eaver endif C Output consts write(*,*)'rlam1=',rlam1,',rlam2=',rlam2,',rlam3=',rlam3 write(10,*)'rlam1=',rlam1,',rlam2=',rlam2,',rlam3=',rlam3 C---------------------------------------------------------------------------- endif !hot and cold start C... Iteration C... mnstep=50 do 210 nstep=istep,mnstep write(10,*)'Iteration=',nstep,' ******************************' write(*,*)'Iteration=',nstep,' ******************************' if(ibifur.eq.1.and.nstep.eq.1) then miter=20 else miter=1 endif do 208 iter=1,miter write(*,*)'rlam1=',rlam1 call coarsen(icsflag,rlam2,rlam3,nd_de,ne_de,nsnp) write(10,*)'Node,elem. deleted=',nd_de,ne_de,',Snapflag=',nsnp write(10,*)'# of nodes=',nn write(*,*)'Node,elem. deleted=',nd_de,ne_de,',Snapflag=',nsnp write(*,*)'# of nodes=',nn call mvnode(1,moved) write(10,*)'Nodes moved=',moved write(*,*)'Nodes moved=',moved call refine(nstep,ibifur,icsflag,rlam1,nd_ex,ne_ex,iover) C Automatic rlam1-3 if(ibifur.eq.1.and.nstep.eq.1) then if(iover.eq.0) then write(*,*)'# of nodes=',nn write(*,*)'Target interval=[',ndmin,ndmax,']' else write(*,*)'# of nodes >',mnp endif if(iter.eq.1) then if(iover.eq.1.or.nn.gt.ndmax) then iul=1 !upper range search rlam1_1=rlam1 rlam1_2=rat_u*rlam1 rlam1=rlam1_2 rlam2=rlam1/4 rlam3=0.8d0*rlam2 call readin(idum) go to 208 else if(nn.lt.ndmin) then iul=0 !lower range search rlam1_1=rat_l*rlam1 rlam1_2=rlam1 rlam1=rlam1_1 rlam2=rlam1/4 rlam3=0.8d0*rlam2 call readin(idum) go to 208 endif endif !iter.eq.1 C 2nd trial if(iter.eq.2) then if(iul.eq.1.and.(iover.eq.1.or.nn.gt.ndmax)) then write(*,*)'Upper automatic search failed' stop endif if(iul.eq.0.and.nn.lt.ndmin) then write(*,*)'Lower automatic search failed' stop endif endif if(iover.eq.0.and.nn.ge.ndmin.and.nn.le.ndmax) then write(*,*)'rlam1-3 found!' write(10,*)'rlam1-3 found!' write(*,*)'rlam1=',rlam1,',rlam2=',rlam2,',rlam3=',rlam3 write(10,*)'rlam1=',rlam1,',rlam2=',rlam2,',rlam3=',rlam3 go to 223 endif C write(*,*)'iover=',iover,',nn=',nn,',ndmin=',ndmin C Bi-section if(iover.eq.1.or.nn.gt.ndmax) then rlam1_1=rlam1 rlam1_2=rlam1_2 else if(nn.lt.ndmin) then rlam1_1=rlam1_1 rlam1_2=rlam1 endif rlam1=(rlam1_1+rlam1_2)/2 rlam2=rlam1/4 rlam3=0.8d0*rlam2 call readin(idum) endif !ibifur.eq.1.and.nstep.eq.1 208 continue !iter=1,miter if(ibifur.eq.1.and.nstep.eq.1) then write(*,*)'Autosearch didnot converge in ',miter,'iterations' stop endif 223 continue write(10,*)'Node/elem. added=',nd_ex,ne_ex write(10,*)'# of nodes=',nn write(*,*)'Node/elem. added=',nd_ex,ne_ex write(*,*)'# of nodes=',nn call mvnode(1,moved) write(10,*)'Nodes moved=',moved write(*,*)'Nodes moved=',moved call swap(nswapped,sperc) write(10,*)'Edge swapped=',nswapped write(*,*)'Edge swapped=',nswapped call mvnode(1,moved) write(10,*)'Nodes moved=',moved write(*,*)'Nodes moved=',moved call store(nstep+1) if(iconflag.eq.0.and.moved.eq.0.and.nswapped.eq.0.and. & nd_ex.eq.0.and.ne_de.eq.0.and.nsnap.eq.0) then write(10,*)'Error optimization completed!' write(*,*)'Error optimization completed!' go to 212 else if(iconflag.ne.0.and.dfloat(nd_de)/nn.lt.dn_r.and. & dfloat(nd_ex)/nn.lt.an_r.and.sperc.lt.se_r) then write(10,*)'Error optimization completed!' write(*,*)'Error optimization completed!' go to 212 endif 210 continue !nstep=istep,mnstep write(*,*)'Failed to converge in ',mnstep,' iterations' stop 212 continue C... C... Output optimized grid C... open(20,file='op.gr3') write(20,*)'Optimized grid in x,y' write(20,*)nt,nn do i=1,nn write(20,100)i,xn(i),yn(i),0 enddo do i=1,nt write(20,101)i,3,(ic1(i,j),j=1,3) enddo close(20) call vscan call integrity write(*,*)'done integrity check after optimization' write(10,*)'done integrity check after optimization' * ********************************************************************************* * * * Geometry run to clean up skew elements (optional) * * Global arrays itab, nrsv are no longer updated; * * a new global array iskew is added * * * ********************************************************************************* * * if(iflaggeo.eq.1) then write(10,*)'Starting geometric run' write(*,*)'Starting geometric run' C Const to determine negative area cst_neg=cst_neg/100 C Initialize skew flags do i=1,nt if(ar(ic1(i,1),ic1(i,2),ic1(i,3)).gt.askew) iskew(i)=1 enddo do it=1,it_geo write(10,*)'Iteration in clean-up=',it,'*********************' write(*,*)'Iteration in clean-up=',it,'*********************' call coarsen_geo(askew,nd_de,nsnap) write(*,*)'Nodes deleted=',nd_de,',Snap flag=',nsnap write(10,*)'Nodes deleted=',nd_de,',Snap flag=',nsnap call vscan call integrity write(*,*)'done integrity check' write(10,*)'done integrity check' call mvnode_geo(askew,10,moved) write(*,*)'Nodes moved=',moved write(10,*)'Nodes moved=',moved call swap_geo(askew,nswapped) write(*,*)'Edge swapped=',nswapped write(10,*)'Edge swapped=',nswapped call mvnode_geo(askew,10,moved) write(*,*)'Nodes moved=',moved write(10,*)'Nodes moved=',moved C Count skew elements nskew=0 do i=1,nt if(iskew(i).eq.1) nskew=nskew+1 enddo !i write(*,*)nskew,' skew elem. left' write(10,*)nskew,' skew elem. left' if(nskew.eq.0) go to 211 if(nd_de.eq.0.and.nsnap.eq.0.and.moved.eq.0.and.nswapped.eq.0) & go to 215 enddo !it=1,it_geo 215 continue write(*,*)'There are still',nskew, ' skew elem. left' write(*,*)'Check fort.10' write(10,*)'There are still',nskew, ' skew elem. left' nskew=0 do i=1,nt if(iskew(i).eq.1) then nskew=nskew+1 write(10,*)nskew,i endif enddo 211 continue endif !iflaggeo.eq.1 * ***************************************************************** * * * End clean-up * * * ***************************************************************** * C... C... Output final grid & reserved element list C... open(20,file='final.gr3') open(25,file='rsv_elem.out') write(20,*)'Final grid in x,y' write(20,*)nt,nn do i=1,nn write(20,100)i,xn(i),yn(i),0 enddo do i=1,nt write(20,101)i,3,(ic1(i,j),j=1,3) enddo close(20) if(ics.eq.2) then C Output in lat and long open(20,file='final.ll') write(20,*)'Final grid in lat-long' write(20,*)nt,nn do i=1,nn C Convert to lat and long. coordinates xll=xn(i)/rearth/dcos(sfea0)+slam0 yll=yn(i)/rearth xll=xll/pi*180 yll=yll/pi*180 write(20,*)i,xll,yll,0 enddo do i=1,nt write(20,101)i,3,(ic1(i,j),j=1,3) enddo close(20) endif !ics=2 if(ires_elem.eq.1) then do i=1,nt write(25,*)i,irsve(i) enddo endif call vscan call integrity write(*,*)'done final integrity check' write(10,*)'done final integrity check' write(10,*)'Completed!' write(*,*)'Completed!' 100 format(i7,2(1x,e24.16),1x,f9.3) 101 format(i6,1x,i1,1x,3i7) stop end * ********************************************************************************* * * * Subroutines * * * ********************************************************************************* * * ********************************************************************************* * * * Store grid for hot start or auto search * * * ********************************************************************************* * subroutine store(nstep) include 'agrid.h' open(20,file='store.dat') rewind(20) write(20,*)nt,nn,nstep do i=1,nn write(20,99)xn(i),yn(i) enddo write(20,*)(nrsv(i),ibdf(i),itab(i),nsur_e(i), & (ic_e(i,j),j=1,nsur_e(i)),i=1,nn) write(20,*)(irsve(i),(ic1(i,j),j=1,3),i=1,nt) close(20) 99 format(2(1x,e34.27)) return end * ********************************************************************************* * * * Read in last stored grid * * * ********************************************************************************* * subroutine readin(istep) include 'agrid.h' open(unit=20,file='store.dat',status='old') read(20,*)nt,nn,istep do i=1,nn read(20,99)xn(i),yn(i) enddo read(20,*)(nrsv(i),ibdf(i),itab(i),nsur_e(i), & (ic_e(i,j),j=1,nsur_e(i)),i=1,nn) read(20,*)(irsve(i),(ic1(i,j),j=1,3),i=1,nt) close(20) 99 format(2(1x,e34.27)) return end * ********************************************************************************* * * * Output intemediate grid * * * ********************************************************************************* * subroutine output_grid include 'agrid.h' open(unit=16,file='tmp.gr3') rewind(16) write(16,*)'Intermediate grid' write(16,*)nt,nn do i=1,nn write(16,*)i,xn(i),yn(i),0 enddo do i=1,nt write(16,*)i,3,(ic1(i,j),j=1,3) enddo close(16) return end * ********************************************************************************* * * * Check validity of the edit grid * * * ********************************************************************************* * subroutine vscan include 'agrid.h' C Check the connectivity do i=1,nt do j=1,3 n1=ic1(i,j) n2=ic1(i,mod(j,3)+1) icount=0 do k=1,nsur_e(n1) ke=ic_e(n1,k) if(ke.ne.i.and. &(ic1(ke,1).eq.n2.or.ic1(ke,2).eq.n2.or.ic1(ke,3).eq.n2)) then icount=icount+1 if(icount.gt.1) then write(*,*)'Error in virus scan; following nodes' write(*,*)'are shared by more than 2 elements:' write(*,*)nd1,nd2 stop endif endif enddo !k=1,nsur_e(n1) if(icount.eq.0.and.(ibdf(n1).eq.0.or.ibdf(n1).eq.0)) then write(*,*)'Fake bnd edge in vscan',n1,n2 call output_grid stop endif enddo !j=1,3 enddo !i=1,nt return end * ********************************************************************************* * * * Integrity check for bnd of edit grid (more comprehensive than vscan) * * * ********************************************************************************* * subroutine integrity include 'agrid.h' dimension isur(mnei),ibnode(2),ibside(mnp,2),iwork(2) C... Find all bnd edges nedg=0 do i=1,nn if(ibdf(i).eq.1) then call pball(i,nsur,isur,nbnode,ibnode) if(nbnode.lt.2.or. & nbnode.eq.2.and.ibnode(1).eq.ibnode(2)) then write(*,*)'Ill-formed bnd at',i write(*,*)nbnode,(ibnode(j),j=1,nbnode) stop endif C Normal situation (nbnode=2) C Check redundancy do 61 k=1,2 iwork(k)=0 !flag do j=1,nedg if(i.eq.ibside(j,1).and.ibnode(k).eq.ibside(j,2).or. & i.eq.ibside(j,2).and.ibnode(k).eq.ibside(j,1)) then iwork(k)=1 go to 61 endif enddo !j 61 continue !k=1,2 do k=1,2 if(iwork(k).eq.0) then nedg=nedg+1 if(nedg.gt.mnp) then write(*,*)'More bnd edges than mnp' stop endif ibside(nedg,1)=i ibside(nedg,2)=ibnode(k) endif enddo !k=1,2 endif !ibdf(i).eq.1 enddo !i=1,nn do i=1,nedg n1=ibside(i,1) n2=ibside(i,2) do 226 j=i+1,nedg n3=ibside(j,1) n4=ibside(j,2) if((n1-n3)*(n1-n4).eq.0.or.(n2-n3)*(n2-n4).eq.0) go to 226 call intersect2(xn(n1),xn(n2),xn(n3),xn(n4), & yn(n1),yn(n2),yn(n3),yn(n4),ifl,xin,yin,tt1,tt2) if(ifl.eq.1) then write(*,*)'Looped boundary at nodes',n1,n2,n3,n4 write(10,*)'Looped boundary at nodes',n1,n2,n3,n4 endif 226 continue !j=i+1,nedg enddo !i=1,nedg return end * ********************************************************************************* * * * Output edge error distribution * * * ********************************************************************************* * subroutine err_distri2(aver,sd,err_min,err_max) include 'agrid.h' dimension eer(3*mne) C Find all edges sum=0 nedg=0 do i=1,nt do j=1,3 nedg=nedg+1 nd1=ic1(i,j) nd2=ic1(i,mod(j,3)+1) rl=error(nd1,nd2,ter) eer(nedg)=ter sum=sum+ter if(nedg.eq.1) then err_min=ter err_max=ter else err_max=dmax1(err_max,ter) err_min=dmin1(err_min,ter) endif enddo !j=1,3 enddo !i=1,nt aver=sum/nedg C Standard deviation open(20,file='err_dis.dat') rewind(20) sum=0 do i=1,nedg sum=sum+(eer(i)-aver)**2 write(20,*)i,eer(i) enddo sd=dsqrt(sum/nedg) close(20) return end * ****************************************************************************** * * * Function to find error for edge (nd1,nd2) (in the edit grid) * * * ****************************************************************************** * function error(nd1,nd2,ter) include 'agrid.h' parameter(nfen=4*2) dimension rint(0:nfen) rl=dsqrt((xn(nd1)-xn(nd2))**2+(yn(nd1)-yn(nd2))**2) if(rl.eq.0) then write(*,*)'Zero edge length at',x1,y1,x2,y2 call output_grid stop endif indicator=0 do i=-1,nfen-1 if(i.le.0) then !compute 2 end nodes first xx=xn(nd1)*(-i)+xn(nd2)*(i+1) yy=yn(nd1)*(-i)+yn(nd2)*(i+1) nnel=itab(nd1)*(-i)+itab(nd2)*(i+1) rint((i+1)*nfen)=err_int(nnel,xx,yy) else xx=xn(nd1)+(xn(nd2)-xn(nd1))/nfen*i yy=yn(nd1)+(yn(nd2)-yn(nd1))/nfen*i if(i.le.nfen/2) then call qsearch(itab(nd1),xx,yy,nnel,nfl) else call qsearch(itab(nd2),xx,yy,nnel,nfl) endif if(nfl.eq.1) then indicator=1 go to 400 endif rint(i)=err_int(nnel,xx,yy) endif enddo !i=-1,nfen-1 400 continue if(indicator.eq.1) then ter=(rint(0)+rint(nfen))/2*rl else ter=dnumingl(rint,nfen,rl) endif error=ter/rl return end * ************************************************************************* * * * Move nodes * * moved: # of nodes moved during last sweeping * * * ************************************************************************* * subroutine mvnode(iter,moved) include 'agrid.h' parameter(theta=0.5) !relaxation const. dimension nnei(mnp),inei(mnp,mnei),erri2(mnei) dimension irsve_p(mnp) C... Surrounding nodes of a node (excluding itself) C... do i=1,nn nnei(i)=0 irsve_p(i)=0 !not in reserved elements do j=1,nsur_e(i) ie=ic_e(i,j) if(irsve(ie).eq.1) irsve_p(i)=1 do 10 k=1,3 nd=ic1(ie,k) if(nd.eq.i) go to 10 C Check redundancy do l=1,nnei(i) if(nd.eq.inei(i,l)) go to 10 enddo !l nnei(i)=nnei(i)+1 if(nnei(i).gt.mnei) then write(*,*)'Too many neigh. in move' stop endif inei(i,nnei(i))=nd 10 continue !k=1,3 enddo !j=1,nsur_e(i) enddo !i=1,nn C... Iteration to move nodes C... do it=1,iter moved=0 do 202 i=1,nn C Check if bnd node (including reserved node), or inside a reserved element if(ibdf(i).eq.1.or.irsve_p(i).eq.1) go to 202 C Store original values xtmp=xn(i) ytmp=yn(i) aerr=0 avx=0 avy=0 do k=1,nnei(i) ndk=inei(i,k) rik=dsqrt((xn(i)-xn(ndk))**2+(yn(i)-yn(ndk))**2) if(rik.eq.0) then write(*,*)'Node i has been moved to its neighbor',i write(*,*)(inei(i,kk),kk=1,nnei(i)) call output_grid stop endif erri2(k)=error(i,ndk,ter)**2 aerr=aerr+erri2(k)*rik*rik/nnei(i) avx=avx+erri2(k)*(xn(i)-xn(ndk))/nnei(i) avy=avy+erri2(k)*(yn(i)-yn(ndk))/nnei(i) enddo !k ss1=0 ss2=0 ss3=0 ss4=0 ss5=0 ss6=0 do j=1,nnei(i) ndj=inei(i,j) rij=dsqrt((xn(i)-xn(ndj))**2+(yn(i)-yn(ndj))**2) rlam1=2*erri2(j)*(xn(i)-xn(ndj))-2*avx rlam2=2*erri2(j)*(yn(i)-yn(ndj))-2*avy ss1=ss1+rlam1**2 ss2=ss2+rlam1*rlam2 ss3=ss3-rlam1*(erri2(j)*rij**2-aerr) ss4=ss2 ss5=ss5+rlam2**2 ss6=ss6-rlam2*(erri2(j)*rij**2-aerr) enddo !j del=ss1*ss5-ss2*ss4 if(del.eq.0) go to 202 del1=ss3*ss5-ss2*ss6 del2=ss1*ss6-ss3*ss4 xn(i)=xn(i)+del1/del*theta yn(i)=yn(i)+del2/del*theta C Check negative area do k=1,nsur_e(i) ne=ic_e(i,k) n1=ic1(ne,1) n2=ic1(ne,2) n3=ic1(ne,3) if(area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)). & lt.cst_neg) then xn(i)=xtmp yn(i)=ytmp go to 202 endif enddo !k moved=moved+1 C Update search table call qsearch(itab(i),xn(i),yn(i),nnel,nfl) if(nfl.eq.1) then itab(i)=itab(i) else itab(i)=nnel endif 202 continue !i=1,nn if(moved.eq.0) go to 204 enddo !it=1,iter 204 continue return end * ************************************************************************* * * * Swap edges * * * ************************************************************************* * subroutine swap(nswap,sperc) include 'agrid.h' dimension asp(2),iwork(2) nswap=0 do 203 i=1,nt if(irsve(i).eq.1) go to 203 do 204 j=1,3 j_1=mod(j,3)+1 j_2=6-j-j_1 node1=ic1(i,j) node2=ic1(i,j_1) call elem2(node1,node2,ie1,ie2) iwork(1)=i iwork(2)=ie1+ie2-i if(iwork(2).eq.0.or.irsve(iwork(2)).eq.1) go to 204 C Compute A.R. do k=1,2 ne=iwork(k) nd1=ic1(ne,1) nd2=ic1(ne,2) nd3=ic1(ne,3) a1=error(nd1,nd2,ter1) a2=error(nd2,nd3,ter2) a3=error(nd3,nd1,ter3) a1=ter1 a2=ter2 a3=ter3 ss=(a1+a2+a3)/2 if(ss.le.a1.or.ss.le.a2.or.ss.le.a3) go to 204 asp(k)=a1*a2*a3/8/(ss-a1)/(ss-a2)/(ss-a3) enddo !k=1,2 C Re-arrange iwork(2) ne=iwork(2) ic1(ne,j_2)=(ic1(ne,1)+ic1(ne,2)+ic1(ne,3))-node1-node2 ic1(ne,j)=node2 ic1(ne,j_1)=node1 C Compute aspect ratios for the 2 alternative elements n11=ic1(iwork(1),j_2) n12=ic1(iwork(1),j) n13=ic1(iwork(2),j_2) n21=ic1(iwork(2),j_2) n22=ic1(iwork(2),j) n23=ic1(iwork(1),j_2) C Check negative area area1=area(xn(n11),xn(n12),xn(n13),yn(n11),yn(n12),yn(n13)) area2=area(xn(n21),xn(n22),xn(n23),yn(n21),yn(n22),yn(n23)) if(area1.lt.cst_neg.or.area2.lt.cst_neg) go to 204 a11=error(n11,n12,ter1) a12=error(n12,n13,ter2) a13=error(n13,n11,ter3) a11=ter1 a12=ter2 a13=ter3 ss1=(a11+a12+a13)/2 if(ss1.le.a11.or.ss1.le.a12.or.ss1.le.a13) go to 204 asp1=a11*a12*a13/8/(ss1-a11)/(ss1-a12)/(ss1-a13) a21=error(n21,n22,ter1) a22=error(n22,n23,ter2) a23=error(n23,n21,ter3) a21=ter1 a22=ter2 a23=ter3 ss2=(a21+a22+a23)/2 if(ss2.le.a21.or.ss2.le.a22.or.ss2.le.a23) go to 204 asp2=a21*a22*a23/8/(ss2-a21)/(ss2-a22)/(ss2-a23) C Compare and swap c if(asp1+asp2.lt.asp(1)+asp(2)) then if(asp1.lt.asp(1).and.asp2.lt.asp(2).or.asp1.lt.asp(2). & and.asp2.lt.asp(1)) then nswap=nswap+1 ic1(iwork(1),j)=n13 ic1(iwork(1),j_1)=n11 ic1(iwork(1),j_2)=n12 ic1(iwork(2),j)=n23 ic1(iwork(2),j_1)=n21 ic1(iwork(2),j_2)=n22 C Update ic_e; only the 4 vertices of the diamond are affected if(nsur_e(n11)+1.gt.mnei.or.nsur_e(n21)+1.gt.mnei) then write(*,*)'Too many neigh. in swap',n11,n21 call output_grid stop endif nsur_e(n11)=nsur_e(n11)+1 ic_e(n11,nsur_e(n11))=iwork(2) nsur_e(n21)=nsur_e(n21)+1 ic_e(n21,nsur_e(n21))=iwork(1) index=0 do l=1,nsur_e(n12) if(ic_e(n12,l).eq.iwork(2)) then index=1 ic_e(n12,l)=ic_e(n12,nsur_e(n12)) go to 18 endif enddo 18 if(index.eq.0) then write(*,*)'Wrong neigh. in swap 1' stop endif nsur_e(n12)=nsur_e(n12)-1 index=0 do l=1,nsur_e(n22) if(ic_e(n22,l).eq.i) then index=1 ic_e(n22,l)=ic_e(n22,nsur_e(n22)) go to 19 endif enddo 19 if(index.eq.0) then write(*,*)'Wrong neigh. in swap 2' stop endif nsur_e(n22)=nsur_e(n22)-1 C Some edges may be left out due to changes in the local indices; C but this can be rectified in the next iteration. endif !A.R. are smaller 204 continue !j=1,3 203 continue !i=1,nt sperc=nswap/(1.5*nt) !estimate return end * ****************************************************************************** * * * Refine edges * * * ****************************************************************************** * subroutine refine(nstep,ibifur,nflag,thres,nd_ex,ne_ex,iover) include 'agrid.h' iover=0 !flag for auto search overflow C... Define theshold if(nflag.eq.0) then C Find the average error sum=0 do i=1,nt do j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) rl=error(k1,k2,ter) sum=sum+ter enddo enddo rlm=sum/3/nt criter=thres*rlm else criter=thres endif C... Mesh refinement nd_ex=0 ne_ex=0 iter=0 i=0 2 continue i=i+1 !element # if(i.gt.nt) go to 5 iter=iter+1 if(iter.gt.2*mne) then write(*,*)'Maximum # of iterations exceeded in refinement!' stop endif if(irsve(i).eq.1) go to 2 C Check extremely small elements n1=ic1(i,1) n2=ic1(i,2) n3=ic1(i,3) if(area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3))/2. & lt.cst_neg) go to 2 do 3 j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) rl=error(k1,k2,ter) if(ter.gt.criter) then C Search for the other element with the side (k1,k2) call elem2(k1,k2,ie1,ie2) i1=ie1+ie2-i if(i1.ne.0.and.irsve(i1).eq.1) go to 3 C Check extremely small elements i1 if(i1.ne.0) then n1=ic1(i1,1) n2=ic1(i1,2) n3=ic1(i1,3) if(area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3))/2. & lt.cst_neg) go to 3 endif if(i1.eq.0) then C Refine bnd edge k3=ic1(i,1)+ic1(i,2)+ic1(i,3)-k1-k2 if(nsur_e(k3)+1.gt.mnei) then write(*,*)'Too many neighbors during refinement' write(*,*)'For auto search,' write(*,*)'try to close up the gap between emin and emax' write(*,*)'Or use manual search instead' write(*,*)'Trouble edge=',k1,k2 call output_grid stop endif C Check overflow if(ibifur.eq.1.and.nstep.eq.1) then if(nn+1.gt.mnp.or.nt+1.gt.mne) then iover=1 go to 5 endif endif if(nn+1.gt.mnp.or.nt+1.gt.mne) then write(*,*)'Too many added nodes/elements in ref.!' call output_grid stop endif nn=nn+1 nd_ex=nd_ex+1 xn(nn)=(xn(k1)+xn(k2))/2 yn(nn)=(yn(k1)+yn(k2))/2 ibdf(nn)=1 nrsv(nn)=0 call qsearch(itab(k1),xn(nn),yn(nn),nnel,nfl) if(nfl.eq.1) then itab(nn)=itab(k1) else itab(nn)=nnel endif nt=nt+1 ne_ex=ne_ex+1 j_1=mod(j,3)+1 ic1(i,j_1)=nn ic1(nt,j)=nn ic1(nt,j_1)=k2 ic1(nt,6-j-j_1)=k3 irsve(nt)=0 C Update ic_e; only k2, k3 and nn are affected nsur_e(nn)=2 ic_e(nn,1)=i ic_e(nn,2)=nt nsur_e(k3)=nsur_e(k3)+1 !bounds checked ic_e(k3,nsur_e(k3))=nt index=0 do l=1,nsur_e(k2) if(ic_e(k2,l).eq.i) then index=1 ic_e(k2,l)=nt go to 9 endif enddo !l 9 if(index.eq.0) then write(*,*)'Wrong neigh. in refine',k2 call output_grid stop endif else !i1/=0 C Refine internal edge k3=ic1(i,1)+ic1(i,2)+ic1(i,3)-k1-k2 k4=ic1(i1,1)+ic1(i1,2)+ic1(i1,3)-k1-k2 if(nsur_e(k3)+1.gt.mnei.or.nsur_e(k4)+1.gt.mnei) then write(*,*)'Too many neighb. in refining intern. edge' write(*,*)'For auto search,' write(*,*)'try to close up the gap between emin and emax' write(*,*)'Or use manual search instead' write(*,*)'Trouble edge=',k1,k2 call output_grid stop endif if(ibifur.eq.1.and.nstep.eq.1) then if(nn+1.gt.mnp.or.nt+2.gt.mne) then iover=1 go to 5 endif endif if(nn+1.gt.mnp.or.nt+2.gt.mne) then write(*,*)'Too many nodes/elements in ref. intern. edge' call output_grid stop endif nn=nn+1 nd_ex=nd_ex+1 xn(nn)=(xn(k1)+xn(k2))/2 yn(nn)=(yn(k1)+yn(k2))/2 ibdf(nn)=0 nrsv(nn)=0 call qsearch(itab(k1),xn(nn),yn(nn),nnel,nfl) if(nfl.eq.1) then itab(nn)=itab(k1) else itab(nn)=nnel endif nt=nt+2 ne_ex=ne_ex+2 j_1=mod(j,3)+1 j_2=6-j-j_1 ic1(i,j_1)=nn ic1(nt-1,j)=nn ic1(nt-1,j_1)=k2 ic1(nt-1,j_2)=k3 irsve(nt-1)=0 ic1(i1,j)=nn ic1(i1,j_1)=k1 ic1(i1,j_2)=k4 ic1(nt,j)=k2 ic1(nt,j_1)=nn ic1(nt,j_2)=k4 irsve(nt)=0 C Update ic_e; only k2-k4 and nn are affected nsur_e(nn)=4 ic_e(nn,1)=i ic_e(nn,2)=i1 ic_e(nn,3)=nt ic_e(nn,4)=nt-1 nsur_e(k3)=nsur_e(k3)+1 ic_e(k3,nsur_e(k3))=nt-1 nsur_e(k4)=nsur_e(k4)+1 ic_e(k4,nsur_e(k4))=nt index=0 do l=1,nsur_e(k2) if(ic_e(k2,l).eq.i) then index=index+1 ic_e(k2,l)=nt-1 endif if(ic_e(k2,l).eq.i1) then index=index+1 ic_e(k2,l)=nt endif enddo !l if(index.ne.2) then write(*,*)'Wrong neigh. in refine internal',k2 call output_grid stop endif endif !i1 C Recheck this element for the next loop i=i-1 go to 2 endif !ter.gt.criter 3 continue !j=1,3 go to 2 5 continue return end * ************************************************************************* * * * Coarsen and snap * * Use of "red" nodes/elements: mark them in the process and delete * * them altogether at the end of coarsen or snap. * * All other connectivity info must be correct and updated for * * all non-red nodes/elements. * * Method for debugging: go thru all with the problem checked * * along the way. * * * ************************************************************************* * subroutine coarsen(nflag,thres1,thres2,nd_de,ne_de,nsnap) include 'agrid.h' dimension iwork(0:mnp),iredp(mnp),irede(mne),nred(0:mne+1) dimension isur(mnei),ibnode1(2),ibnode2(2) *-----------------------------------------------------------------------* * * * Coarsen edges * * * *-----------------------------------------------------------------------* do i=1,nn iredp(i)=0 !initialize enddo do i=1,nt irede(i)=0 enddo C... Define theshold if(nflag.eq.0) then C Find the average error sum=0 do i=1,nt do j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) rl=error(k1,k2,ter) sum=sum+ter enddo enddo rlm=sum/3/nt criter=thres1*rlm else criter=thres1 endif C... Nodes merging C... do 21 i=1,nt if(irede(i).eq.1.or.irsve(i).eq.1) go to 21 do 31 j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) rl=error(k1,k2,ter) if(ter.lt.criter) then C Search for the other element with the side (k1,k2) call elem2(k1,k2,ie1,ie2) i1=ie1+ie2-i if(i1.ne.0.and.irsve(i1).eq.1) go to 31 C... Check if both are reserved nodes if(nrsv(k1).eq.1.and.nrsv(k2).eq.1) go to 31 C... Check if it's a narrow region if(ibdf(k1).eq.1.and.ibdf(k2).eq.1.and.i1.ne.0) then c write(11,*)'Narrow region warning!' go to 31 endif C... Check if collapsing bnd (k1,k2) would lead to topological change; C... this happends when they are among the 2 of 3-node bnd if(i1.eq.0) then call pball(k1,nsur,isur,nbnode1,ibnode1) call pball(k2,nsur,isur,nbnode2,ibnode2) if(nbnode1.ne.2.or.nbnode2.ne.2) then write(*,*)'Bad bnd ball in csn',nbnode1,nbnode2 stop endif k7=ibnode1(1)+ibnode1(2)-k2 k8=ibnode2(1)+ibnode2(2)-k1 if(ibdf(k7).ne.1.or.ibdf(k8).ne.1) then write(*,*)'Wrong bnd ball in csn' call output_grid stop endif if(k7.eq.k8) go to 31 endif C... Define k3, k4 C... if(ibdf(k1).ne.ibdf(k2)) then k4=k1*ibdf(k1)+k2*ibdf(k2) k3=(k1+k2)-k4 else if(nrsv(k1).eq.1.or.nrsv(k2).eq.1) then k4=k1*nrsv(k1)+k2*nrsv(k2) k3=(k1+k2)-k4 else !k2->k1 k3=k2 k4=k1 endif C... Move k3 --> k4 C... Check negative areas and if reserved C... do 12 l=1,nsur_e(k3) nel=ic_e(k3,l) if(irsve(nel).eq.1) go to 31 if(nel.eq.i.or.nel.eq.i1) go to 12 do lin=1,3 if(ic1(nel,lin).eq.k3) then l_1=mod(lin,3)+1 n1=k4 n2=ic1(nel,l_1) n3=ic1(nel,6-lin-l_1) if(area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)). & lt.cst_neg) go to 31 go to 12 endif enddo !lin=1,3 12 continue !l do 15 l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.eq.i.or.nel.eq.i1) go to 15 do lin=1,3 if(ic1(nel,lin).eq.k3) ic1(nel,lin)=k4 enddo 15 continue C... Stack red nodes/elements C... iredp(k3)=1 irede(i)=1 if(i1.ne.0) irede(i1)=1 C... Update neighborhood C... Only the 4 vertices of the diamond are affected. do lin=1,3 nd=ic1(i,lin) index=0 do l=1,nsur_e(nd) if(ic_e(nd,l).eq.i) then index=1 ic_e(nd,l)=ic_e(nd,nsur_e(nd)) go to 16 endif enddo 16 if(index.eq.0) then write(*,*)'Wrong neigh. in csn',nd,k3,k4,i call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !lin=1,3 if(i1.ne.0) then do lin=1,3 nd=ic1(i1,lin) index=0 do l=1,nsur_e(nd) if(ic_e(nd,l).eq.i1) then index=1 ic_e(nd,l)=ic_e(nd,nsur_e(nd)) go to 17 endif enddo 17 if(index.eq.0) then write(*,*)'Wrong neigh. csn 2',nd,k3,k4,i1 call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !lin=1,3 endif !i1.ne.0 C Lump k3's ball to k4's do l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.eq.i.or.nel.eq.i1) then write(*,*)'Ball of k3 has been updated',i,i1 stop endif nsur_e(k4)=nsur_e(k4)+1 if(nsur_e(k4).gt.mnei) then write(*,*)'Too many neighbors in csn',k4 stop endif ic_e(k4,nsur_e(k4))=nel enddo nsur_e(k3)=0 go to 21 !element i has been deleted endif !rl.lt.criter 31 continue !j=1,3 21 continue !i=1,nt C... Delete all red elements and update reserved element list C... C Order all red elements nred_m=0 !# do i=1,nt if(irede(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nt+1 !for convenience do i=1,nred_m do j=nred(i)+1,nred(i+1)-1 do k=1,3 ic1(j-i,k)=ic1(j,k) enddo !k irsve(j-i)=irsve(j) enddo !j enddo !i nt=nt-nred_m ne_de=nred_m C... Delete hanging nodes and update bnd flags, reserved node flags, and search table C... nred_m=0 do i=1,nn if(iredp(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nn+1 do i=1,nred_m do k=nred(i)+1,nred(i+1)-1 xn(k-i)=xn(k) yn(k-i)=yn(k) nrsv(k-i)=nrsv(k) itab(k-i)=itab(k) C Bnd nodes may have been deleted; but internal and bnd nodes do not change C their nature. ibdf(k-i)=ibdf(k) enddo enddo C Generate a search table to speed up next section nred(0)=0 !For convenience do i=0,nred_m do k=nred(i)+1,nred(i+1)-1 iwork(k)=i enddo if(i.ne.0) iwork(nred(i))=-1 !Flag enddo do l=1,nt do lin=1,3 if(iwork(ic1(l,lin)).eq.-1) then write(*,*)'Fake hanging node',ic1(l,lin) call output_grid stop endif ic1(l,lin)=ic1(l,lin)-iwork(ic1(l,lin)) enddo enddo nn=nn-nred_m nd_de=nred_m C... Update neighborhood call surrounding(1) *-----------------------------------------------------------------------* * * * Snap nodes * * * *-----------------------------------------------------------------------* C Const to indicate if this stage has been performed; used to C prevent the loop in main programme been terminated while C this stage is still on nsnap=0 C... No red nodes in this module do i=1,nt irede(i)=0 enddo C... Define theshold if(nflag.eq.0) then sum=0 do i=1,nt do j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) rl=error(k1,k2,ter) sum=sum+ter enddo enddo rlm=sum/3/nt criter=thres2*rlm else criter=thres2 endif C... Snap C... do 22 i=1,nt if(irsve(i).eq.1) go to 22 dism=-100 do j=1,3 j_1=mod(j,3)+1 l1=ic1(i,j) l2=ic1(i,j_1) rl=dsqrt((xn(l1)-xn(l2))**2+(yn(l1)-yn(l2))**2) if(rl.gt.dism) then dism=rl k1=l1 k2=l2 k3=ic1(i,6-j-j_1) j0=j j01=j_1 endif enddo !j=1,3 C... Estimate snapping "distance" xx=(xn(k1)+xn(k2))/2 yy=(yn(k1)+yn(k2))/2 dist_mid=dsqrt((xx-xn(k3))**2+(yy-yn(k3))**2) err_mid=(err_int(itab(k1),xn(k1),yn(k1))+ + err_int(itab(k2),xn(k2),yn(k2)))/2 rl=(err_int(itab(k3),xn(k3),yn(k3))+err_mid)/2*dist_mid if(rl.lt.criter) then C Find neighboring elements of (k1,k2) call elem2(k1,k2,ie1,ie2) i2=ie1+ie2-i if(i2.ne.0.and.irsve(i2).eq.1) go to 22 if(i2.eq.0) then C Bnd edge; shift k3 to (xx,yy) if(ibdf(k3).eq.1) go to 22 C Check negative areas and if reserved do l=1,nsur_e(k3) nel=ic_e(k3,l) if(irsve(nel).eq.1) go to 22 if(nel.ne.i) then do lin=1,3 if(ic1(nel,lin).eq.k3) then l_1=mod(lin,3)+1 n2=ic1(nel,l_1) n3=ic1(nel,6-lin-l_1) if(area(xx,xn(n2),xn(n3),yy,yn(n2),yn(n3)). & lt.cst_neg) go to 22 endif enddo !lin endif !nel.ne.i enddo !l=1,nsur_e(k3) C Update search table call qsearch(itab(k3),xx,yy,nnel,nfl) if(nfl.eq.1) then itab(k3)=itab(k3) else itab(k3)=nnel endif xn(k3)=xx yn(k3)=yy C Stack red element irede(i)=1 C Update ic_e; only k1-k3 are affected do l=1,3 nd=ic1(i,l) !element i not eliminated yet index=0 do ll=1,nsur_e(nd) if(ic_e(nd,ll).eq.i) then index=1 ic_e(nd,ll)=ic_e(nd,nsur_e(nd)) go to 11 endif enddo !ll 11 if(index.eq.0) then write(*,*)'Wrong neigh. in snap',nd,i call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !l=1,3 C One more bnd node ibdf(k3)=1 C Update the status const nsnap=1 else !i2/=0 C Internal edge; swap edge C Re-arrange the local indices of 2 elements iwork(1)=i iwork(2)=i2 do l=1,2 ne=iwork(l) nd1=ic1(ne,1) nd2=ic1(ne,2) nd3=ic1(ne,3) npro=isign(1,nd1-nd2)*isign(1,nd2-nd3)*isign(1,nd3-nd1) nd1p=(nd1+nd2+nd3)-(k1+k2) nd2p=k1 nd3p=k2 npro3=isign(1,nd1p-nd2p)*isign(1,nd2p-nd3p)* *isign(1,nd3p-nd1p) if(npro3.eq.npro) then ic1(ne,1)=nd1p ic1(ne,2)=nd2p ic1(ne,3)=nd3p else if(npro3.eq.-npro) then ic1(ne,1)=nd1p ic1(ne,2)=nd3p ic1(ne,3)=nd2p else write(*,*)'Error' stop endif enddo !l=1,2 C 2 alternative elements n11=ic1(iwork(1),1) n12=ic1(iwork(1),2) n13=ic1(iwork(2),1) n21=ic1(iwork(2),1) n22=ic1(iwork(2),2) n23=ic1(iwork(1),1) C Check negative area a1=area(xn(n11),xn(n12),xn(n13),yn(n11),yn(n12),yn(n13)) a2=area(xn(n21),xn(n22),xn(n23),yn(n21),yn(n22),yn(n23)) if(a1.lt.cst_neg.or.a2.lt.cst_neg) go to 22 C Swap ic1(iwork(1),1)=n11 ic1(iwork(1),2)=n12 ic1(iwork(1),3)=n13 ic1(iwork(2),1)=n21 ic1(iwork(2),2)=n22 ic1(iwork(2),3)=n23 C Update ic_e; only the 4 vertices of the diamond are affected if(nsur_e(n11)+1.gt.mnei.or.nsur_e(n21)+1.gt.mnei) then write(*,*)'Too many neigh. in snap',n11,n21 call output_grid stop endif nsur_e(n11)=nsur_e(n11)+1 ic_e(n11,nsur_e(n11))=i2 nsur_e(n21)=nsur_e(n21)+1 ic_e(n21,nsur_e(n21))=i index=0 do l=1,nsur_e(n12) if(ic_e(n12,l).eq.i2) then index=1 ic_e(n12,l)=ic_e(n12,nsur_e(n12)) go to 18 endif enddo 18 if(index.eq.0) then write(*,*)'Wrong neigh. in snap 1' call output_grid stop endif nsur_e(n12)=nsur_e(n12)-1 index=0 do l=1,nsur_e(n22) if(ic_e(n22,l).eq.i) then index=1 ic_e(n22,l)=ic_e(n22,nsur_e(n22)) go to 19 endif enddo 19 if(index.eq.0) then write(*,*)'Wrong neigh. in snap 2' call output_grid stop endif nsur_e(n22)=nsur_e(n22)-1 C Update the status const nsnap=1 endif !i2=0 endif !rl.lt.criter 22 continue !i=1,nt C... Delete all red elements nred_m=0 !# do i=1,nt if(irede(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nt+1 !for convenience do i=1,nred_m do j=nred(i)+1,nred(i+1)-1 do k=1,3 ic1(j-i,k)=ic1(j,k) enddo !k irsve(j-i)=irsve(j) enddo !j enddo !i nt=nt-nred_m ne_de=ne_de+nred_m C.. Update neighborhood call surrounding(1) return end * ************************************************************************* * * * Clean-up module 1: Move nodes * * moved: # of nodes moved during last sweeping * * * ************************************************************************* * subroutine mvnode_geo(askew,iter,moved) include 'agrid.h' parameter(theta=0.5) !relaxation const. dimension nnei(mnp),inei(mnp,mnei),irsve_p(mnp),asp(mnei) C... Surrounding nodes of a node (excluding itself) C... do i=1,nn nnei(i)=0 irsve_p(i)=0 !not in reserved elements do j=1,nsur_e(i) ie=ic_e(i,j) if(irsve(ie).eq.1) irsve_p(i)=1 do 10 k=1,3 nd=ic1(ie,k) if(nd.eq.i) go to 10 C Check redundancy do l=1,nnei(i) if(nd.eq.inei(i,l)) go to 10 enddo !l nnei(i)=nnei(i)+1 if(nnei(i).gt.mnei) then write(*,*)'Too many neigh. in geo_move' stop endif inei(i,nnei(i))=nd 10 continue !k=1,3 enddo !j=1,nsur_e(i) enddo !i=1,nn C... Iteration to improve A.R. C... do it=1,iter moved=0 do 202 ielem=1,nt if(iskew(ielem).eq.0.or.irsve(ielem).eq.1) go to 202 do 207 jj=1,3 i=ic1(ielem,jj) C Check if bnd node (including reserved node) or inside reserved elem. if(ibdf(i).eq.1.or.irsve_p(i).eq.1) go to 207 C Store original values and compute original A.R. xtmp=xn(i) ytmp=yn(i) do k=1,nsur_e(i) ne=ic_e(i,k) asp(k)=ar(ic1(ne,1),ic1(ne,2),ic1(ne,3)) enddo !k aerr=0 avx=0 avy=0 do k=1,nnei(i) ndk=inei(i,k) rik=dsqrt((xn(i)-xn(ndk))**2+(yn(i)-yn(ndk))**2) if(rik.eq.0) then write(*,*)'Node i has been moved to its neigh. in geo',i call output_grid stop endif aerr=aerr+rik*rik/nnei(i) avx=avx+(xn(i)-xn(ndk))/nnei(i) avy=avy+(yn(i)-yn(ndk))/nnei(i) enddo !k ss1=0 ss2=0 ss3=0 ss4=0 ss5=0 ss6=0 do j=1,nnei(i) ndj=inei(i,j) rij=dsqrt((xn(i)-xn(ndj))**2+(yn(i)-yn(ndj))**2) rlam1=2*(xn(i)-xn(ndj))-2*avx rlam2=2*(yn(i)-yn(ndj))-2*avy ss1=ss1+rlam1**2 ss2=ss2+rlam1*rlam2 ss3=ss3-rlam1*(rij**2-aerr) ss4=ss2 ss5=ss5+rlam2**2 ss6=ss6-rlam2*(rij**2-aerr) enddo !j del=ss1*ss5-ss2*ss4 if(del.eq.0) go to 207 del1=ss3*ss5-ss2*ss6 del2=ss1*ss6-ss3*ss4 xn(i)=xn(i)+del1/del*theta yn(i)=yn(i)+del2/del*theta C Check negative area and A.R. do k=1,nsur_e(i) ne=ic_e(i,k) n1=ic1(ne,1) n2=ic1(ne,2) n3=ic1(ne,3) if(area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)). & lt.cst_neg.or. & ar(n1,n2,n3).gt.asp(k).and.ar(n1,n2,n3).gt.askew) then xn(i)=xtmp yn(i)=ytmp go to 207 endif enddo !k moved=moved+1 C Update skew element list do k=1,nsur_e(i) ne=ic_e(i,k) if(ar(ic1(ne,1),ic1(ne,2),ic1(ne,3)).le.askew) then iskew(ne)=0 else iskew(ne)=1 endif enddo !k 207 continue !jj=1,3 202 continue !ielem=1,nt if(moved.eq.0) go to 204 enddo !it=1,iter 204 continue return end * ************************************************************************* * * * Clean-up module 2: Swap edges * * * ************************************************************************* * subroutine swap_geo(askew,nswap) include 'agrid.h' dimension asp(2),iwork(2) nswap=0 do 203 i=1,nt if(irsve(i).eq.1.or.iskew(i).eq.0) go to 203 do 204 j=1,3 j_1=mod(j,3)+1 j_2=6-j-j_1 node1=ic1(i,j) node2=ic1(i,j_1) call elem2(node1,node2,ie1,ie2) iwork(1)=i iwork(2)=ie1+ie2-i if(iwork(2).eq.0.or.irsve(iwork(2)).eq.1) go to 204 C Compute A.R. do k=1,2 ne=iwork(k) asp(k)=ar(ic1(ne,1),ic1(ne,2),ic1(ne,3)) enddo !k=1,2 C Re-arrange iwork(2) ne=iwork(2) ic1(ne,j_2)=(ic1(ne,1)+ic1(ne,2)+ic1(ne,3))-node1-node2 ic1(ne,j)=node2 ic1(ne,j_1)=node1 C Compute aspect ratios for the 2 alternative elements n11=ic1(iwork(1),j_2) n12=ic1(iwork(1),j) n13=ic1(iwork(2),j_2) n21=ic1(iwork(2),j_2) n22=ic1(iwork(2),j) n23=ic1(iwork(1),j_2) C Check negative area area1=area(xn(n11),xn(n12),xn(n13),yn(n11),yn(n12),yn(n13)) area2=area(xn(n21),xn(n22),xn(n23),yn(n21),yn(n22),yn(n23)) if(area1.lt.cst_neg.or.area2.lt.cst_neg) go to 204 asp1=ar(n11,n12,n13) asp2=ar(n21,n22,n23) C Compare and swap if(asp1.le.askew.and.asp2.le.askew.or. & asp1.lt.asp(1).and.asp2.lt.asp(2).or. & asp1.lt.asp(2).and.asp2.lt.asp(1)) then nswap=nswap+1 ic1(iwork(1),j)=n13 ic1(iwork(1),j_1)=n11 ic1(iwork(1),j_2)=n12 ic1(iwork(2),j)=n23 ic1(iwork(2),j_1)=n21 ic1(iwork(2),j_2)=n22 C Update ic_e; only the 4 vertices of the diamond are affected if(nsur_e(n11)+1.gt.mnei.or.nsur_e(n21)+1.gt.mnei) then write(*,*)'Too many neigh. in swap',n11,n21 call output_grid stop endif nsur_e(n11)=nsur_e(n11)+1 ic_e(n11,nsur_e(n11))=iwork(2) nsur_e(n21)=nsur_e(n21)+1 ic_e(n21,nsur_e(n21))=iwork(1) index=0 do l=1,nsur_e(n12) if(ic_e(n12,l).eq.iwork(2)) then index=1 ic_e(n12,l)=ic_e(n12,nsur_e(n12)) go to 18 endif enddo 18 if(index.eq.0) then write(*,*)'Wrong neigh. in geo_swap 1' stop endif nsur_e(n12)=nsur_e(n12)-1 index=0 do l=1,nsur_e(n22) if(ic_e(n22,l).eq.i) then index=1 ic_e(n22,l)=ic_e(n22,nsur_e(n22)) go to 19 endif enddo 19 if(index.eq.0) then write(*,*)'Wrong neigh. in geo_swap 2' stop endif nsur_e(n22)=nsur_e(n22)-1 C Update skew element list if(asp1.le.askew) then iskew(iwork(1))=0 else iskew(iwork(1))=1 endif if(asp2.le.askew) then iskew(iwork(2))=0 else iskew(iwork(2))=1 endif endif !A.R. are smaller if(asp1.le.askew.and.asp2.le.askew) go to 203 204 continue !j=1,3 203 continue !i=1,nt return end * ************************************************************************* * * * Clean-up module 3: Coarsen and snap * * * ************************************************************************* * subroutine coarsen_geo(askew,nd_de,nsnap) include 'agrid.h' dimension iwork(0:mnp),iredp(mnp),irede(mne),nred(0:mne+1) dimension isur(mnei),ibnode1(2),ibnode2(2) *-----------------------------------------------------------------------* * * * Take out 3-node balls * * * *-----------------------------------------------------------------------* do i=1,nn iredp(i)=0 !initialize enddo do i=1,nt irede(i)=0 enddo C... Nodes merging C... i=0 21 continue i=i+1 !elem. # if(i.gt.nt) go to 601 if(irede(i).eq.1.or.irsve(i).eq.1.or.iskew(i).eq.0) go to 21 C See if it can be removed by combining with other 2 elements to C form a single triangle do 51 j=1,3 n1=ic1(i,j) if(ibdf(n1).eq.0.and.nsur_e(n1).eq.3) then j_1=mod(j,3)+1 n2=ic1(i,j_1) n3=ic1(i,6-j-j_1) call elem2(n1,n2,ie1,ie2) i2=ie1+ie2-i call elem2(n1,n3,ie1,ie2) i3=ie1+ie2-i if(i2.eq.0.or.i3.eq.0) then write(*,*)'Inconsis. ball in geo_csn' stop endif if(irsve(i2).eq.1.or.irsve(i3).eq.1) go to 51 n4=ic1(i2,1)+ic1(i2,2)+ic1(i2,3)-n1-n2 if(area(xn(n2),xn(n3),xn(n4),yn(n2),yn(n3),yn(n4)). & lt.cst_neg) then write(*,*)'Wrong orie. in 3-ball node' call output_grid stop endif if(ar(n2,n3,n4).le.ar(n1,n2,n3)) then C Take out n1 nt=nt+1 if(nt.gt.mne) then write(*,*)'Too many elements in geo_csn' stop endif ic1(nt,1)=n2 ic1(nt,2)=n3 ic1(nt,3)=n4 irsve(nt)=0 if(ar(n2,n3,n4).le.askew) then iskew(nt)=0 else iskew(nt)=1 endif iredp(n1)=1 iwork(1)=i iwork(2)=i2 iwork(3)=i3 do k=1,3 ie=iwork(k) irede(ie)=1 iskew(ie)=0 !unimportant C Neighborhood do kk=1,3 nd=ic1(ie,kk) index=0 do l=1,nsur_e(nd) if(ic_e(nd,l).eq.ie) then index=1 ic_e(nd,l)=ic_e(nd,nsur_e(nd)) go to 52 endif enddo !l 52 if(index.eq.0) then write(*,*)'Wrong neigh. in geo_csn',nd,ie call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !kk=1,3 enddo !k=1,3 iwork(1)=n2 iwork(2)=n3 iwork(3)=n4 do k=1,3 nd=iwork(k) nsur_e(nd)=nsur_e(nd)+1 !bound checked if(nsur_e(nd).gt.mnei) then write(*,*)'Impossible in geo_csn',nd stop endif ic_e(nd,nsur_e(nd))=nt enddo !k go to 21 endif !ar(n2,n3,n4).le.ar(n1,n2,n3) endif !ibdf(n1).eq.0.and.nsur_e(n1).eq.3 51 continue !j=1,3 601 continue C... C... Delete all red elements and update reserved and skew element lists C... C Order all red elements nred_m=0 !# do i=1,nt if(irede(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nt+1 !for convenience do i=1,nred_m do j=nred(i)+1,nred(i+1)-1 do k=1,3 ic1(j-i,k)=ic1(j,k) enddo !k irsve(j-i)=irsve(j) iskew(j-i)=iskew(j) enddo !j enddo !i nt=nt-nred_m C... Delete hanging nodes and update bnd flags C... nred_m=0 do i=1,nn if(iredp(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nn+1 do i=1,nred_m do k=nred(i)+1,nred(i+1)-1 xn(k-i)=xn(k) yn(k-i)=yn(k) C Bnd nodes may have been deleted; but internal and bnd nodes do not change C their nature. ibdf(k-i)=ibdf(k) enddo enddo C Generate a search table to speed up next section nred(0)=0 !For convenience do i=0,nred_m do k=nred(i)+1,nred(i+1)-1 iwork(k)=i enddo if(i.ne.0) iwork(nred(i))=-1 !Flag enddo do l=1,nt do lin=1,3 if(iwork(ic1(l,lin)).eq.-1) then write(*,*)'Fake hanging node in geo',ic1(l,lin) call output_grid stop endif ic1(l,lin)=ic1(l,lin)-iwork(ic1(l,lin)) enddo enddo nn=nn-nred_m nd_de=nred_m C... Update neighborhood call surrounding(1) call vscan call integrity write(*,*)'Passed virus scan after extracting 3-node balls' *-----------------------------------------------------------------------* * * * Coarsen edges * * * *-----------------------------------------------------------------------* C... do i=1,nn iredp(i)=0 !initialize enddo do i=1,nt irede(i)=0 enddo C... Nodes merging C... i=0 23 continue i=i+1 !elem. # if(i.gt.nt) go to 602 if(irede(i).eq.1.or.irsve(i).eq.1.or.iskew(i).eq.0) go to 23 C Look for an edge to coarsen do 31 j=1,3 k1=ic1(i,j) k2=ic1(i,mod(j,3)+1) C Search for the other element with the side (k1,k2) call elem2(k1,k2,ie1,ie2) i1=ie1+ie2-i if(i1.ne.0.and.irsve(i1).eq.1) go to 31 C... Check if both are reserved nodes c if(nrsv(k1).eq.1.and.nrsv(k2).eq.1) go to 31 C... Check if it's a narrow region if(ibdf(k1).eq.1.and.ibdf(k2).eq.1.and.i1.ne.0) then c write(11,*)'Narrow region warning!' go to 31 endif C... Check if collapsing bnd (k1,k2) would lead to topological change; C... this happends when they are among the 2 of 3-node bnd if(i1.eq.0) then call pball(k1,nsur,isur,nbnode1,ibnode1) call pball(k2,nsur,isur,nbnode2,ibnode2) if(nbnode1.ne.2.or.nbnode2.ne.2) then write(*,*)'Bad bnd edge in geo_csn',nbnode1,nbnode2 call output_grid stop endif k7=ibnode1(1)+ibnode1(2)-k2 k8=ibnode2(1)+ibnode2(2)-k1 if(ibdf(k7).ne.1.or.ibdf(k8).ne.1) then write(*,*)'Wrong bnd ball in geo_csn' call output_grid stop endif if(k7.eq.k8) go to 31 endif !i1.eq.0 C... Shift nodes k1, k2 to their new position k4 C... if(ibdf(k1).eq.ibdf(k2)) then !k2->k1 k3=k2 k4=k1 else k4=k1*ibdf(k1)+k2*ibdf(k2) k3=(k1+k2)-k4 endif C... Move k3 --> k4 C... Check negative areas, A.R. and if reserved C... Swap k3 and k4 if unsuccessful (but k1, k2 must be of same type) C... icount=0 42 continue do 12 l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.eq.i.or.nel.eq.i1) go to 12 asp0=ar(ic1(nel,1),ic1(nel,2),ic1(nel,3)) do lin=1,3 if(ic1(nel,lin).eq.k3) then l_1=mod(lin,3)+1 n1=k4 n2=ic1(nel,l_1) n3=ic1(nel,6-lin-l_1) areal=area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)) if(irsve(nel).eq.1.or. & areal.lt.cst_neg.or.areal.ge.cst_neg.and. & asp0.le.askew.and.ar(n1,n2,n3).gt.askew) then c & ar(n1,n2,n3).gt.asp0.and.ar(n1,n2,n3).gt.askew) then if(icount.eq.0.and.ibdf(k1).eq.ibdf(k2)) then C Swap k3, k4 and retry ktmp=k3 k3=k4 k4=ktmp icount=icount+1 go to 42 else go to 31 endif endif !irsve(nel).eq.1.or.... go to 12 endif !ic1(nel,lin).eq.k3 enddo !lin=1,3 12 continue !l=1,nsur_e(k3) do 15 l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.eq.i.or.nel.eq.i1) go to 15 do lin=1,3 if(ic1(nel,lin).eq.k3) ic1(nel,lin)=k4 enddo 15 continue C... Stack red nodes/elements C... iredp(k3)=1 irede(i)=1 if(i1.ne.0) irede(i1)=1 C... Update neighborhood C... Only the 4 vertices of the diamond are affected. C... do lin=1,3 nd=ic1(i,lin) index=0 do l=1,nsur_e(nd) if(ic_e(nd,l).eq.i) then index=1 ic_e(nd,l)=ic_e(nd,nsur_e(nd)) go to 16 endif enddo 16 if(index.eq.0) then write(*,*)'Wrong neigh. in geo_csn',nd,k3,k4,i call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !lin=1,3 if(i1.ne.0) then do lin=1,3 nd=ic1(i1,lin) index=0 do l=1,nsur_e(nd) if(ic_e(nd,l).eq.i1) then index=1 ic_e(nd,l)=ic_e(nd,nsur_e(nd)) go to 17 endif enddo 17 if(index.eq.0) then write(*,*)'Wrong neigh. geo_csn 2',nd,k3,k4,i1 call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !lin=1,3 endif !i1.ne.0 C Lump k3's ball to k4's do l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.eq.i.or.nel.eq.i1) then write(*,*)'Ball of k3 has been updated 2',i,i1 stop endif nsur_e(k4)=nsur_e(k4)+1 if(nsur_e(k4).gt.mnei) then write(*,*)'Too many neighbors in geo_csn',k4 stop endif ic_e(k4,nsur_e(k4))=nel enddo nsur_e(k3)=0 C Update iskew iskew(i)=0 !unimportant if(i1.ne.0) iskew(i1)=0 do l=1,nsur_e(k4) nel=ic_e(k4,l) if(ar(ic1(nel,1),ic1(nel,2),ic1(nel,3)).gt.askew) then iskew(nel)=1 else iskew(nel)=0 endif enddo !l go to 23 !element i has been deleted 31 continue !j=1,3 go to 23 602 continue C... C... Delete all red elements and update reserved and skew element lists C... C Order all red elements nred_m=0 !# do i=1,nt if(irede(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nt+1 !for convenience do i=1,nred_m do j=nred(i)+1,nred(i+1)-1 do k=1,3 ic1(j-i,k)=ic1(j,k) enddo !k irsve(j-i)=irsve(j) iskew(j-i)=iskew(j) enddo !j enddo !i nt=nt-nred_m C... Delete hanging nodes and update bnd flags C... nred_m=0 do i=1,nn if(iredp(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nn+1 do i=1,nred_m do k=nred(i)+1,nred(i+1)-1 xn(k-i)=xn(k) yn(k-i)=yn(k) C Bnd nodes may have been deleted; but internal and bnd nodes do not change C their nature. ibdf(k-i)=ibdf(k) enddo enddo C Generate a search table to speed up next section nred(0)=0 !For convenience do i=0,nred_m do k=nred(i)+1,nred(i+1)-1 iwork(k)=i enddo if(i.ne.0) iwork(nred(i))=-1 !Flag enddo do l=1,nt do lin=1,3 if(iwork(ic1(l,lin)).eq.-1) then write(*,*)'Fake hanging node in geo',ic1(l,lin) call output_grid stop endif ic1(l,lin)=ic1(l,lin)-iwork(ic1(l,lin)) enddo enddo nn=nn-nred_m nd_de=nd_de+nred_m C... Update neighborhood call surrounding(1) call vscan call integrity write(*,*)'Passed virus scan after coarsening' *-----------------------------------------------------------------------* * * * Snap nodes to bnd edges (swap internal edge redundant with swap_geo) * * * *-----------------------------------------------------------------------* C Const to indicate if this stage has been performed; used to C prevent the loop in main programme been terminated while C this stage is still on nsnap=0 C... No red nodes in this module do i=1,nt irede(i)=0 enddo C... Snap C... do 22 i=1,nt if(irsve(i).eq.1.or.iskew(i).eq.0) go to 22 C Try to snap to eliminate i dism=-100 do j=1,3 j_1=mod(j,3)+1 l1=ic1(i,j) l2=ic1(i,j_1) rl=dsqrt((xn(l1)-xn(l2))**2+(yn(l1)-yn(l2))**2) if(rl.gt.dism) then dism=rl k1=l1 k2=l2 k3=ic1(i,6-j-j_1) j0=j j01=j_1 endif enddo !j=1,3 C Foot of projection if(dism.eq.0) then write(*,*)'Degenerate element in geo_csn',i stop endif tt=((xn(k3)-xn(k1))*(xn(k2)-xn(k1))+ +(yn(k3)-yn(k1))*(yn(k2)-yn(k1)))/dism/dism if(tt.le.0.or.tt.ge.1) then write(*,*)'Foot outside of (k1,k2) in geo_csn' stop endif xx=xn(k1)+(xn(k2)-xn(k1))*tt yy=yn(k1)+(yn(k2)-yn(k1))*tt C Find neighboring elements of (k1,k2) call elem2(k1,k2,ie1,ie2) i2=ie1+ie2-i if(i2.ne.0) go to 22 C Snap bnd edge; shift k3 to (xx,yy) if(ibdf(k3).eq.1) go to 22 C Check negative areas, A.R. and if reserved do l=1,nsur_e(k3) nel=ic_e(k3,l) if(nel.ne.i) then asp0=ar(ic1(nel,1),ic1(nel,2),ic1(nel,3)) do lin=1,3 if(ic1(nel,lin).eq.k3) then l_1=mod(lin,3)+1 n2=ic1(nel,l_1) n3=ic1(nel,6-lin-l_1) areal=area(xx,xn(n2),xn(n3),yy,yn(n2),yn(n3)) a1=dsqrt((xx-xn(n2))**2+(yy-yn(n2))**2) a2=dsqrt((xn(n2)-xn(n3))**2+(yn(n2)-yn(n3))**2) a3=dsqrt((xn(n3)-xx)**2+(yn(n3)-yy)**2) if(areal.gt.0) asp=dmax1(a1,a2,a3)/dsqrt(areal/pi) if(irsve(nel).eq.1.or.areal.lt.cst_neg.or. & areal.gt.cst_neg.and.asp0.le.askew.and.asp.gt.askew) go to 22 endif enddo !lin=1,3 endif !nel.ne.i enddo !l=1,nsur_e(k3) xn(k3)=xx yn(k3)=yy C Stack red element irede(i)=1 C Update ic_e; only k1-k3 are affected do l=1,3 nd=ic1(i,l) !element i not eliminated yet index=0 do ll=1,nsur_e(nd) if(ic_e(nd,ll).eq.i) then index=1 ic_e(nd,ll)=ic_e(nd,nsur_e(nd)) go to 11 endif enddo !ll 11 if(index.eq.0) then write(*,*)'Wrong neigh. in geo_snap',nd,i call output_grid stop endif nsur_e(nd)=nsur_e(nd)-1 enddo !l=1,3 C One more bnd node ibdf(k3)=1 C Update iskew iskew(i)=0 !unimportant do l=1,nsur_e(k3) nel=ic_e(k3,l) if(ar(ic1(nel,1),ic1(nel,2),ic1(nel,3)).gt.askew) then iskew(nel)=1 else iskew(nel)=0 endif enddo !l C Update the status const nsnap=1 22 continue !i=1,nt C... Delete all red elements nred_m=0 !# do i=1,nt if(irede(i).eq.1) then nred_m=nred_m+1 nred(nred_m)=i endif enddo nred(nred_m+1)=nt+1 !for convenience do i=1,nred_m do j=nred(i)+1,nred(i+1)-1 do k=1,3 ic1(j-i,k)=ic1(j,k) enddo !k irsve(j-i)=irsve(j) iskew(j-i)=iskew(j) enddo !j enddo !i nt=nt-nred_m C.. Update neighborhood call surrounding(1) call vscan call integrity write(*,*)'Passed virus scan after swapping' return end * ********************************************************************************* * * * Compute ball of a node * * igrid=0: bg grid; =1: edit grid. * * ic_e is distinctive from each other. * * * ********************************************************************************* * subroutine surrounding(igrid) include 'agrid.h' if(igrid.eq.0) then np=nn0 ne=nt0 else np=nn ne=nt endif do i=1,np nsur_e(i)=0 enddo do i=1,ne do j=1,3 if(igrid.eq.0) then nd=ic0(i,j) else nd=ic1(i,j) endif nsur_e(nd)=nsur_e(nd)+1 if(nsur_e(nd).gt.mnei) then if(igrid.eq.0) then write(*,*)'Too many neighbors in bg grid' else write(*,*)'Too many neighbors in edit grid' call output_grid endif stop endif ic_e(nd,nsur_e(nd))=i enddo enddo return end * ********************************************************************************* * * * Create 2-tier neighborhood of a node in bg grid. * * ic3 includes the node itself * * * ********************************************************************************* * subroutine surrounding2 include 'agrid.h' dimension iwork((mnei+1)*mnei) call surrounding(0) do i=1,nn0 nsur2(i)=1 ic3(i,1)=i C Find all elements do j=1,nsur_e(i) iwork(j)=ic_e(i,j) enddo icount=nsur_e(i) do j=1,nsur_e(i) do 12 k=1,3 ndc=ic0(ic_e(i,j),k) if(ndc.eq.i) go to 12 do 10 l=1,nsur_e(ndc) ie=ic_e(ndc,l) C Check redundancy do m=1,icount if(ie.eq.iwork(m)) go to 10 enddo !m icount=icount+1 !bound guaranteed if(icount.gt.(mnei+1)*mnei) then write(*,*)'Impossible in 2-tier ball' stop endif iwork(icount)=ie 10 continue !l=1,nsur_e(ndc) 12 continue !k=1,3 enddo !j=1,nsur_e(i) C Find all nodes do j=1,icount do 11 k=1,3 nd=ic0(iwork(j),k) C Check redundancy do l=1,nsur2(i) if(nd.eq.ic3(i,l)) go to 11 enddo !l nsur2(i)=nsur2(i)+1 if(nsur2(i).gt.mnei) then write(*,*)'Too many neigh. in 2-tier ball' stop endif ic3(i,nsur2(i))=nd 11 continue !k=1,3 enddo !j=1,icount enddo !i=1,nn0 return end * ********************************************************************************* * * * Least-square filter for error field. * * * ********************************************************************************* * subroutine efilter(error1) include 'agrid.h' dimension error1(mnp),coe(6,7),rhs(6) do 10 i=1,nn0 if(nsur2(i).le.6) then error1(i)=error0(i) go to 10 endif s40=0 s41=0 s42=0 s43=0 s44=0 s30=0 s31=0 s32=0 s33=0 s20=0 s21=0 s22=0 s10=0 s11=0 s00=0 rhs(1)=0 rhs(2)=0 rhs(3)=0 rhs(4)=0 rhs(5)=0 rhs(6)=0 do j=1,nsur2(i) nd=ic3(i,j) !including i xlc=x0(nd) ylc=y0(nd) C Coefficients s40=s40+xlc**4 !x^4*y^0 s41=s41+xlc**3*ylc s42=s42+xlc**2*ylc**2 s43=s43+xlc*ylc**3 s44=s44+ylc**4 s30=s30+xlc**3 s31=s31+xlc**2*ylc s32=s32+xlc*ylc**2 s33=s33+ylc**3 s20=s20+xlc**2 s21=s21+xlc*ylc s22=s22+ylc**2 s10=s10+xlc s11=s11+ylc s00=s00+1 rhs(1)=rhs(1)+xlc**2*error0(nd) rhs(2)=rhs(2)+xlc*ylc*error0(nd) rhs(3)=rhs(3)+ylc**2*error0(nd) rhs(4)=rhs(4)+xlc*error0(nd) rhs(5)=rhs(5)+ylc*error0(nd) rhs(6)=rhs(6)+error0(nd) enddo !j C Augmented matrix coe(1,1)=s40 coe(1,2)=s41 coe(1,3)=s42 coe(1,4)=s30 coe(1,5)=s31 coe(1,6)=s20 coe(2,1)=s41 coe(2,2)=s42 coe(2,3)=s43 coe(2,4)=s31 coe(2,5)=s32 coe(2,6)=s21 coe(3,1)=s42 coe(3,2)=s43 coe(3,3)=s44 coe(3,4)=s32 coe(3,5)=s33 coe(3,6)=s22 coe(4,1)=s30 coe(4,2)=s31 coe(4,3)=s32 coe(4,4)=s20 coe(4,5)=s21 coe(4,6)=s10 coe(5,1)=s31 coe(5,2)=s32 coe(5,3)=s33 coe(5,4)=s21 coe(5,5)=s22 coe(5,6)=s11 coe(6,1)=s20 coe(6,2)=s21 coe(6,3)=s22 coe(6,4)=s10 coe(6,5)=s11 coe(6,6)=s00 coe(1,7)=rhs(1) coe(2,7)=rhs(2) coe(3,7)=rhs(3) coe(4,7)=rhs(4) coe(5,7)=rhs(5) coe(6,7)=rhs(6) call dsid(coe,1,6,6,7,det) if(det.eq.0) then error1(i)=error0(i) else error1(i)=coe(1,7)*x0(i)**2+coe(2,7)*x0(i)*y0(i)+ +coe(3,7)*y0(i)**2+coe(4,7)*x0(i)+coe(5,7)*y0(i)+coe(6,7) if(error1(i).lt.0) error1(i)=error0(i) endif 10 continue !i=1,nn0 return end function area(x1,x2,x3,y1,y2,y3) implicit real*8(a-h,o-z) area=0.5*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) return end * ********************************************************************************* * * * Find error at a pt (xx,yy) knowing its parent element nnel * * If (xx,yy) is not in nnel, the results may not be accurate. * * * ********************************************************************************* * function err_int(nnel,xx,yy) include 'agrid.h' m1=ic0(nnel,1) m2=ic0(nnel,2) m3=ic0(nnel,3) s1=dabs(area(xx,x0(m2),x0(m3),yy,y0(m2),y0(m3))) s2=dabs(area(x0(m1),xx,x0(m3),y0(m1),yy,y0(m3))) s3=dabs(area(x0(m1),x0(m2),xx,y0(m1),y0(m2),yy)) err_int=(error0(m1)*s1+error0(m2)*s2+error0(m3)*s3)/area0(nnel) return end * ********************************************************************************* * * * Find 2 adjacent elements (ie1,ie2) for edge (n1,n2) in the edit grid. * * * ********************************************************************************* * subroutine elem2(n1,n2,ie1,ie2) include 'agrid.h' dimension iwork(2) icount=0 do l=1,nsur_e(n1) ne=ic_e(n1,l) if(ic1(ne,1).eq.n2.or.ic1(ne,2).eq.n2.or.ic1(ne,3).eq.n2) then icount=icount+1 if(icount.gt.2) then write(*,*)'More than 2 elements sharing same edge' write(*,*)n1,n2 call output_grid stop endif iwork(icount)=ne endif enddo if(icount.eq.0) then write(*,*)'Hanging edge',n1,n2 call output_grid stop endif ie1=iwork(1) if(icount.eq.1) then ie2=0 !bnd flag else !icount=2 ie2=iwork(2) endif return end * ********************************************************************************* * * * Find ball of nodes (excluding itself) for a node in edit grid; * * for a bnd node, also output 2 adjacent bnd nodes, using the info on * * ball of elements and bnd flags. * * * ********************************************************************************* * subroutine pball(node,nsur,isur,nbnode,ibnode) include 'agrid.h' dimension isur(mnei),ibnode(2) nsur=0 nbnode=0 do j=1,nsur_e(node) ie=ic_e(node,j) do 10 k=1,3 nd=ic1(ie,k) if(nd.eq.node) go to 10 C Check redundancy do l=1,nsur if(nd.eq.isur(l)) go to 10 enddo !l nsur=nsur+1 if(nsur.gt.mnei) then write(*,*)'Too many neigh. in pball' stop endif isur(nsur)=nd if(ibdf(node).eq.1.and.ibdf(nd).eq.1) then call elem2(node,nd,ie1,ie2) if(ie2.eq.0) then nbnode=nbnode+1 if(nbnode.gt.2) then write(*,*)'> 2 adjecent bnd nds in pball' stop endif ibnode(nbnode)=nd endif endif 10 continue !k=1,3 enddo !j=1,nsur_e(i) return end * ********************************************************************************* * * * Compute A.R. in Euclidean space for 3 ordered nodes (n1,n2,n3) * * in edit grid * * * ********************************************************************************* * function ar(n1,n2,n3) include 'agrid.h' a1=dsqrt((xn(n1)-xn(n2))**2+(yn(n1)-yn(n2))**2) a2=dsqrt((xn(n2)-xn(n3))**2+(yn(n2)-yn(n3))**2) a3=dsqrt((xn(n3)-xn(n1))**2+(yn(n3)-yn(n1))**2) area1=area(xn(n1),xn(n2),xn(n3),yn(n1),yn(n2),yn(n3)) if(area1.le.0) then write(*,*)'Negative area in ar()',n1,n2,n3 stop endif ar=dmax1(a1,a2,a3)/dsqrt(area1/pi) return end * ********************************************************************************* * * * Straightline search algorithm. nel is an estimate of element near * * (xt,yt). nfl=1 if bnd is hit. * * * ********************************************************************************* * subroutine qsearch(nel,xt,yt,nnel,nfl) include 'agrid.h' sums(x1,x2,x3,x4,y1,y2,y3,y4)= & dabs((x4-x3)*(y2-y3)+(x2-x3)*(y3-y4))/2+ + dabs((x4-x1)*(y3-y1)-(y4-y1)*(x3-x1))/2+ + dabs((y4-y1)*(x2-x1)-(x4-x1)*(y2-y1))/2 nfl=0 nel00=nel C Starting element nel n1=ic0(nel,1) n2=ic0(nel,2) n3=ic0(nel,3) x1=x0(n1) x2=x0(n2) x3=x0(n3) y1=y0(n1) y2=y0(n2) y3=y0(n3) x4=xt y4=yt aa=sums(x1,x2,x3,x4,y1,y2,y3,y4) ae=dabs(aa-area0(nel))/area0(nel) if(ae.lt.epsilon) then nnel=nel go to 400 endif C (xt,yt) not in nel, and thus (xcg,ycg) and (xt,yt) are distinctive C add asymetric noise to avoid "tangential" cut xcg=(x1+x2+x3)/3 ycg=(y1+y2+y3)/3 xcg=xcg*0.997+x1*(1-0.997) ycg=ycg*0.997+y1*(1-0.997) xcg=xcg*0.998+x2*(1-0.998) ycg=ycg*0.998+y2*(1-0.998) xcg=xcg*0.9991+x3*(1-0.9991) ycg=ycg*0.9991+y3*(1-0.9991) pathl=dsqrt((xt-xcg)**2+(yt-ycg)**2) if(pathl.eq.0) then write(*,*)'Zero path',x0,y0,xt,yt,xcg,ycg stop endif md1=0 !for first intersecting edge only md2=0 it=0 852 continue it=it+1 if(it.gt.3*mne/2) then write(*,*)'Too many searched edges' write(*,*)nel00,xcg,ycg,xt,yt stop endif C Intersecting edge nel_j do 854 j=1,3 jd1=ic0(nel,j) jd2=ic0(nel,mod(j,3)+1) if(jd1.eq.md1.and.jd2.eq.md2.or. &jd2.eq.md1.and.jd1.eq.md2) go to 854 if(x0(jd1).eq.x0(jd2).and.y0(jd1).eq.y0(jd2)) then write(*,*)'Zero line',jd1,jd2,x0(jd1),y0(jd1),x0(jd2),y0(jd2) stop endif call intersect2(xcg,xt,x0(jd1),x0(jd2), & ycg,yt,y0(jd1),y0(jd2),iflag,xin,yin,tt1,tt2) if(iflag.eq.1) then nel_j=j go to 399 endif 854 continue !j=1,3 write(*,*)'Found no intersecting edges' write(*,*)nel00,nel,md1,md2,xt,yt stop 399 md1=ic0(nel,nel_j) md2=ic0(nel,mod(nel_j,3)+1) C Check exit if(ic4(nel,nel_j).eq.0) then nfl=1 go to 400 endif C Check nel's neighbor with edge nel_j nel=ic4(nel,nel_j) !next front element k1=ic0(nel,1) k2=ic0(nel,2) k3=ic0(nel,3) aa=sums(x0(k1),x0(k2),x0(k3),xt,y0(k1),y0(k2),y0(k3),yt) ae=dabs(aa-area0(nel))/area0(nel) if(ae.lt.epsilon) then nnel=nel go to 400 endif go to 852 400 continue return end * ********************************************************************************* * * * Program to detect if two segments (1,2) and (3,4) have common pts * * Assumption: the 4 pts are distinctive. * * The eqs. for the 2 lines are: X=X1+(X2-X1)*tt1 and X=X3+(X4-X3)*tt2. * * Output: iflag: 0: no intersection or colinear; 1: exactly 1 intersection. * * If iflag=1, (xin,yin) is the intersection. * * * ********************************************************************************* * subroutine intersect2(x1,x2,x3,x4,y1,y2,y3,y4,iflag,xin,yin, &tt1,tt2) implicit real*8(a-h,o-z) parameter(zero=0.0) !small positive number or 0 tt1=-1000 tt2=-1000 iflag=0 delta=(x2-x1)*(y3-y4)-(y2-y1)*(x3-x4) delta1=(x3-x1)*(y3-y4)-(y3-y1)*(x3-x4) delta2=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) if(delta.eq.0.0d0) then C------------------------------------------------------- C Check if the two are colinear c do l=1,4 c if(l.eq.1) then c cosi=angle(x1,x2,x1,x3,y1,y2,y1,y3) c rat=rmag(x1,x2,x1,x3,y1,y2,y1,y3) c else if(l.eq.2) then c cosi=angle(x1,x2,x1,x4,y1,y2,y1,y4) c rat=rmag(x1,x2,x1,x4,y1,y2,y1,y4) c else if(l.eq.3) then c cosi=angle(x3,x4,x3,x1,y3,y4,y3,y1) c rat=rmag(x3,x4,x3,x1,y3,y4,y3,y1) c else if(l.eq.4) then c cosi=angle(x3,x4,x3,x2,y3,y4,y3,y2) c rat=rmag(x3,x4,x3,x2,y3,y4,y3,y2) c endif c if(cosi.eq.1.and.rat.ge.1) iflag=-1 c enddo C------------------------------------------------------- else !delta\=0 tt1=delta1/delta tt2=delta2/delta if(tt1.ge.-zero.and.tt1.le.1+zero.and. &tt2.ge.-zero.and.tt2.le.1+zero) then iflag=1 xin=x1+(x2-x1)*tt1 yin=y1+(y2-y1)*tt1 endif C------------------------------------------------------- endif return end * ********************************************************************************* * * * Standard fortrans * * * ********************************************************************************* * SUBROUTINE DSID(SIDW,LSID,NROW,MSID,NSID,SIDET) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION SIDW(NROW,NSID),IFSID(2000),IGSID(2000),ILSID(2000) C SIDET=1.0 SIDET=1.0D00 ISID1=MSID-1 ISIDR=NSID-MSID IF(ISIDR.LT.0) GO TO 16 DO 1 KSID=1,MSID ILSID(KSID)=0 1 IGSID(KSID)=KSID KSID1=1 DO 8 KSID=1,MSID IF(LSID.GE.0) KSID1=KSID+1 C RSID=0.0 RSID=0.0D+0 DO 2 ISID=1,MSID IF(ILSID(ISID).NE.0) GO TO 2 WSID=SIDW(ISID,KSID) XSID=ABS(WSID) IF(RSID.GT.XSID) GO TO 2 RSID=XSID PSID=WSID KFSID=ISID 2 CONTINUE IFSID(KSID)=KFSID ILSID(KFSID)=KFSID SIDET=SIDET*PSID C IF(SIDET.EQ.0.0) GO TO 16 IF(SIDET.EQ.0.0D+0) GO TO 16 C XSID=1.0/PSID XSID=1.0D+0/PSID DO 4 ISID=1,MSID IF(ISID.NE.KFSID) GO TO 3 SIDW(ISID,KSID)=XSID GO TO 4 3 SIDW(ISID,KSID)=-SIDW(ISID,KSID)*XSID 4 CONTINUE IF(KSID1.GT.NSID) GO TO 8 DO 7 JSID=KSID1,NSID IF(JSID.EQ.KSID) GO TO 7 WSID=SIDW(KFSID,JSID) C IF(WSID.EQ.0.0) GO TO 7 IF(WSID.EQ.0.0D+0) GO TO 7 DO 6 ISID=1,MSID IF(ISID.NE.KFSID) GO TO 5 SIDW(ISID,JSID)=WSID*XSID GO TO 6 5 SIDW(ISID,JSID)=SIDW(ISID,JSID)+WSID*SIDW(ISID,KSID) 6 CONTINUE 7 CONTINUE 8 CONTINUE DO 15 KSID=1,ISID1 KFSID=IFSID(KSID) KLSID=ILSID(KFSID) KGSID=IGSID(KSID) IF(KFSID.EQ.KGSID) GO TO 15 IF(LSID) 10,14,9 9 IF(ISIDR) 16,14,12 10 DO 11 ISID=1,MSID WSID=SIDW(ISID,KGSID) SIDW(ISID,KGSID)=SIDW(ISID,KFSID) 11 SIDW(ISID,KFSID)=WSID 12 DO 13 JSID=1,NSID WSID=SIDW(KLSID,JSID) SIDW(KLSID,JSID)=SIDW(KSID,JSID) 13 SIDW(KSID,JSID)=WSID 14 ILSID(KFSID)=KSID ILSID(KGSID)=KLSID IGSID(KLSID)=IGSID(KSID) IGSID(KSID)=KFSID SIDET=-SIDET 15 CONTINUE 16 RETURN END C Program to solve a cyclic tridiagonal system C a,b,c,r are input vectors that are not modified; they are in C the same sense as in tridag. alpha, beta are the lower and upper C corner entries. x is the unknown. subroutine cyclic(a,b,c,alpha,beta,r,x,n) implicit real*8(a-h,o-z) parameter(nmax=100000) dimension a(n),b(n),c(n),r(n),x(n),bb(nmax),u(nmax),z(nmax) dimension gam(nmax) if(n.le.2) then write(*,*)'n too small in cyclic' stop endif if(n.gt.nmax) then write(*,*)'nmax too small in cyclic' stop endif if(b(1).eq.0.) then write(*,*)'re-arrange the matrix' stop endif gamma=-b(1) !avoid underflow bb(1)=b(1)-gamma !set up the new diagonal bb(n)=b(n)-alpha*beta/gamma do i=2,n-1 bb(i)=b(i) enddo call tridag(a,bb,c,r,x,gam,n) u(1)=gamma u(n)=alpha do i=2,n-1 u(i)=0 enddo call tridag(a,bb,c,u,z,gam,n) fact=(x(1)+beta*x(n)/gamma)/(1+z(1)+beta*z(n)/gamma) do i=1,n x(i)=x(i)-fact*z(i) enddo return end function dnumingl(f,n,xlengt) implicit real*8(a-h,o-z) dimension f(0:n),store(3),tore(2) istep=4 icount=0 200 continue nsteps=n/istep h=xlengt/nsteps icount=icount+1 sum=-(f(0)+f(n))/2.0d0 do 600 i=0,n,istep sum=sum+f(i) 600 continue sum=h*sum store(icount)=sum istep=istep/2 if(istep.gt.0) goto 200 do 700 i=1,2 tore(i)=(4.0d0*store(i+1)-store(i))/3.0d0 700 continue dnumingl=(16.0d0*tore(2)-tore(1))/15.0d0 if(xlengt.lt.0.0d0) then dnumingl=-dnumingl endif c write(6,*)'The answer from using the Trapezoidal rules:' c write(6,*) (store(i),i=1,3) c write(6,*)'The answer from using the Simpsons rules:' c write(6,*) (tore(i),i=1,2) c write(6,*)'The best answer is supposed to be:' c write(6,*) dnumingl return end C This program solves a cyclic tridiagonal system. It was copied and adapted C from "Numerical recipes in FORTRAN (pp.43, 68). subroutine tridag(a,b,c,r,u,gam,n) C Solve a tridiagonal system. C a,b,c,r are input vectors and are not modified by this program. C b is the main diagonal, a below and c above. a(1) and c(n) are C not used. r is the r.h.s., gam is a wroking array, C and u is the unknown. implicit real*8(a-h,o-z) ! parameter(nmax=10000) dimension a(n),b(n),c(n),r(n),u(n),gam(n) if(b(1).eq.0.) pause 'tridag: rewrite equations' C u2 should be eliminated bet=b(1) u(1)=r(1)/bet do j=2,n gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.) pause 'tridag failed' u(j)=(r(j)-a(j)*u(j-1))/bet enddo C Backsubstitution do j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) enddo return end