subroutine boundary( i maxpo,isegs,nloop, i lpbnd,idpnt,lpbnds,idseg,hmax,hmin, i numpo,numseg,npbg,rpbg, i icseg,cseg,nrcseg, m nface,lface,vface,coeg,lseg,ibou,ifarea,lquad, m nquad,rquad,npoin,maxseg) * ************************************************************************* * in: * * maxpo maximum number of points * imesh mesh number * isegs segment data * (1,i): name of i:th segment * (2,i): number of points defining segment * (3,i): segment type * nloop number of closed loops of segments " * lpbnd(*) pointer to point names on segments in the loop * idpnt(*) list of points on segments in loop * lpbnds(*) pointer to segment names in the loop * idseg(*) list of segments in the loop * icseg(4,*) list of curved segments * (1,*): segment name * (2,*): startpoint * (3,*): endpoint * (4,*): number of control points * cseg(*) control points and normal vectors for curved segs * nrcseg number of curved segs * hloc meshisze * numpo number of points on all segments (points on * internal boundaries are counted twice) * numseg number of segments on current area (segments on * internal boundaries are counted twice) * * * ************************************************************************* * modified: * nface number of boundary faces * lface(3,nface) 1: node 1 * 2: node 2 (triangulates left of 1-2) * 3: 0 - inactive, 1,2 - active * vface(4,nface) 1: x-normal point 1 * 2: y-normal point 1 * 3: x-normal point 2 * 4: y-normal point 2 * coeg(2,nno) coordinates of points in the local (u,v)-plane * ibou(nno) 1 for boundary point, -1 for internal point * lseg(2,*) info about the triangulation of (global) seg. * 1: 0- seg. not triangulated, >0- number of faces * 2: start position in lface * ifarea(*) boundary faces * (1): number of faces * (2-): boundary faces * lquad(7,nquad) quadtree of points * nquad number of quads * rquad(3,nquad) geometry definition of quad * npoin number of points * maxseg largest segment number ************************************************************************* * implicit double precision (a-h,o-z) dimension idpnt(*),ibou(*),rpbg(npbg,*),coeg(2,*) dimension lface(3,*),isegs(3,*),lpbnds(*),vface(4,*) dimension lquad(7,*),rquad(3,*),icseg(4,*),cseg(nrcseg,4,*) dimension lseg(2,*),ifarea(*),idseg(*),lpbnd(*),vn(2) common /files/ infile,iofile,iscree * ********************************************************* * set up quad for boundary information ********************************************************* * * call dwmatg(icseg,4,nrcseg,0,'i','icseg') * call dwmatg(cseg,4,2,0,'rd','cseg') * call dwmatg(isegs,3,numseg,0,'i','isegs') * call dwmatg(lpbnds,1,nloop,0,'i','lpbnds') * call dwmatg(lpbnd,1,nloop,0,'i','lpbnd') * call dwmatg(idpnt,1,2*numpo,0,'i','idpnt') * write(6,*) 'nloop=',nloop,' numseg=',numseg * maxseg=0 i1=0 nquad=1 do 11 ipoin=1,numpo ip1=iabs(idpnt(ipoin)) ibou(ip1)=1 call addq(lquad,nquad,rquad,ip1,coeg) * 11 continue * npol=1 nfol=nface+1 nmax=0 * * write(*,*) 'numpo=',numpo do 300 iloop=1,nloop lpoins=lpbnds(iloop) lpoint=lpbnd(iloop) * ********************************************************* * loop over all boundary segments, nloop=numbo ********************************************************* * if(nloop.eq.1) then nsegl=numseg else if(iloop.lt.nloop) then lp2=lpbnds(iloop+1) nsegl=lp2-lpoins else nsegl=numseg-lpoins+1 endif endif * ********************************************************* * segment n:o isname ********************************************************* * icount=0 * call dwmatg(idseg(lpoins),1,nsegl,0,'i','idesg') * call dwmatg(idpnt(lpoint),1,nsegl,0,'i','idpnt') * do 200 ij=1,nsegl * iss=idseg(lpoins+ij-1) isname=isegs(1,iss) istart=nface+1 * icount=icount+1 ip1=idpnt(lpoint+icount-1) ip2=idpnt(lpoint+icount) if(ij.eq.nsegl.and.ip2.ge.0) ip2=idpnt(lpoint) if(ij.eq.nsegl.and.ip2.lt.0) ip2=-ip2 * write(*,*) 'ip1=',ip1,' ip2=',ip2 * iseg=isname maxseg=max0(maxseg,iseg) * ********************************************************* * check if already triangulated, ********************************************************* * if(lseg(1,iseg).gt.0) then inr=lseg(1,iseg) ist=lseg(2,iseg) ien=ist+inr-1 * do 101 j=ien,ist,-1 lface(3,j)=2 101 continue icount=icount-1 goto 200 endif * ********************************************************* * not triangulated, check if curved ********************************************************* * call ichecc(nrcseg,icseg,isname,itrue,inamcs) * ********************************************************* * check whether to change the order of points ********************************************************* * if(itrue.eq.1) then nextra=icseg(4,inamcs) if(icseg(2,inamcs).ne.ip1) then call shuffle( $ cseg,nrcseg,inamcs,nextra) icseg(3,inamcs)=icseg(2,inamcs) icseg(2,inamcs)=ip1 endif ipstep=1 endif * xs=coeg(1,ip1) ys=coeg(2,ip1) * ipoin=ip1 told=0.d0 * 110 continue * ********************************************************* * compute local meshsize ********************************************************* * call shepard( i xs,ys,npbg,rpbg(1,1),rpbg(1,2),rpbg(1,3), o hm) hm=dmax1(hm,hmin) hm=dmin1(hm,hmax) xsi=xs ysi=ys ip=ipoin if(itrue.eq.0) call straight(ip,ip1,ip2,coeg,xs,ys,hm) if(itrue.eq.1) call spline(ip,ipstep,ip1,ip2,coeg,told, $ xs,ys,vn,hm,nextra,cseg,nrcseg,inamcs) * ipold=ipoin ipoin=npoin+1 if(ip.eq.ip2) ipoin=ip2 * nface=nface+1 i1=i1+1 ifarea(i1+1)=nface lface(3,nface)=1 lface(1,nface)=ipold lface(2,nface)=ipoin * * Normal vectors * if(itrue.eq.1) then if(ipold.eq.ip1) then vface(1,nface)=cseg(inamcs,1,1) vface(2,nface)=cseg(inamcs,2,1) else vface(1,nface)=vface(3,nface-1) vface(2,nface)=vface(4,nface-1) endif if(ipoin.eq.ip2) then vface(3,nface)=cseg(inamcs,3,1) vface(4,nface)=cseg(inamcs,4,1) else vface(3,nface)=vn(1) vface(4,nface)=vn(2) endif endif * if(ipoin.eq.ip2) goto 100 * npoin=npoin+1 coeg(1,ipoin)=xs coeg(2,ipoin)=ys ibou(ipoin)=1 * call addq(lquad,nquad,rquad,ipoin,coeg) goto 110 100 continue * ********************************************************* * fill segment info ********************************************************* * lseg(1,iseg)=nface-istart+1 lseg(2,iseg)=istart * 200 continue 300 continue ifarea(1)=i1 write(6,*) ' Number of boundary faces:',i1 write(6,*) ' Number of boundary points:',npoin * call dwmatg(coeg,2,npoin,0,'rd','coeg') * call dwmatg(lface,3,nface,0,'i','lface') * call dwmatg(lseg,2,13,0,'i','lseg') * return end subroutine gridc( & nheap,nface,lface,coeg,lheap,rface,lpoin,lfapo,nfapo, x itnode,lquad,npoin,xyquad,nquad,nt,mele,ierr, x ibou,hmax,hmin,iarea,npbg,rpbg) * implicit double precision (a-h,o-z) dimension xyquad(3,*),coeg(2,*),rface(*),itnode(3,*) dimension lface(3,*),lheap(*),lpoin(*),lfapo(3,*),iqp(3600) dimension lquad(7,*),vect(2),lhp(3600),ikno(3) dimension rhp(3600) dimension dirs(2),xyadd(2,5) dimension ibou(*),rpbg(npbg,*) common /files/ infile,iofile,iscree * eps=1.d-9 * mele=180 * write(6,*) 'max element=' * read(5,*) mele * if(mele.eq.0) then * ierr=12 * return * endif 10 continue * ********************************************************* * select face to be deleted ********************************************************* * if(nheap.eq.0) goto 9000 call remh(lheap,rface,nheap,iface) if(lface(3,iface).eq.0) goto 10 11 continue inod1=lface(1,iface) inod2=lface(2,iface) * xs=coeg(1,inod1) xe=coeg(1,inod2) ys=coeg(2,inod1) ye=coeg(2,inod2) * ********************************************************* * compute local meshsize ********************************************************* * xm=(xs+xe)/2.d0 ym=(ys+ye)/2.d0 call shepard( i xm,ym,npbg,rpbg(1,1),rpbg(1,2),rpbg(1,3), o h1) h1=dmin1(h1,hmax) h1=dmax1(h1,hmin) * ********************************** s1=1.d0 vect(1)=1.d0 vect(2)=0.d0 ******************************** dx=xe-xs dy=ye-ys dl=dsqrt(dx*dx+dy*dy) dirs(1)=-dy/dl dirs(2)=dx/dl ikno(1)=inod1 ikno(2)=inod2 xinod1=coeg(1,inod1) yinod1=coeg(2,inod1) xinod2=coeg(1,inod2) yinod2=coeg(2,inod2) * lfa = lface(3,iface) if(lfa.gt.0) lface(3,iface) = lfa - 1 * * check if this boundary is to be triangulated again * if(lface(3,iface).gt.0) then lfa=lface(1,iface) lface(1,iface)=lface(2,iface) lface(2,iface)=lfa rface(iface)=dl call addh(lheap,rface,nheap,iface) endif * xm=(xs+xe)/2.d0 ym=(ys+ye)/2.d0 * * check direction * if((dirs(1)*vect(1)+dirs(2)*vect(2)).lt.0) then vect(1)=-vect(1) vect(2)=-vect(2) endif * vx=xinod2-xinod1 vy=yinod2-yinod1 dl=dsqrt(vx*vx+vy*vy) dirs(1)=-vy/dl dirs(2)=vx/dl * * h1=hloc/sc * if(h1.gt.1.5d0*dl) h1=1.5d0*dl if(h1.lt.0.55d0*dl) h1=0.55d0*dl xm=(xinod1+xinod2)/2.d0 ym=(yinod1+yinod2)/2.d0 * dlrot=dsqrt(h1*h1-0.25d0*dl*dl) xyadd(1,1)=dirs(1)*dlrot+xm xyadd(2,1)=dirs(2)*dlrot+ym xyadd(1,2)=4.d0*dirs(1)*dlrot/5.d0+xm xyadd(2,2)=4.d0*dirs(2)*dlrot/5.d0+ym xyadd(1,3)=3.d0*dirs(1)*dlrot/5.d0+xm xyadd(2,3)=3.d0*dirs(2)*dlrot/5.d0+ym xyadd(1,4)=2.d0*dirs(1)*dlrot/5.d0+xm xyadd(2,4)=2.d0*dirs(2)*dlrot/5.d0+ym xyadd(1,5)=dirs(1)*dlrot/5.d0+xm xyadd(2,5)=dirs(2)*dlrot/5.d0+ym xp=xyadd(1,1) yp=xyadd(2,1) * * find quad * radi=3.d0*h1 radi=2.d0*h1 * call poinqq(lquad,xyquad,xp,yp,radi,iqp,npq,coeg) * * check all points in iquad * i1=0 if(npq.eq.0) goto 35 do 30 ipq=1,npq ip=iqp(ipq) * * check side * xip=coeg(1,ip) yip=coeg(2,ip) xipa=xip yipa=yip dx=xip-xm dy=yip-ym prod=dirs(1)*dx+dirs(2)*dy if(prod.lt.eps) goto 30 dx=xip-xp dy=yip-yp i1=i1+1 lhp(i1)=ip rhp(i1)=dsqrt(dx*dx+dy*dy) 30 continue call sort(lhp,rhp,i1) * * add nodes closer to the segment * 35 continue do 42 iv=1,5 42 lhp(i1+iv)=npoin+iv itot=i1+5 * ********************************************************* * check order of points to be chosen ********************************************************* * call relist(lhp,rhp,xp,yp,coeg,h1,itot) * * choose node * do 50 i=1,itot inod3=lhp(i) * ikno(3)=inod3 * ********************************************************* * check that no face is crossed ********************************************************* * call facech(lhp,npoin,itot,lface,lpoin,lfapo,inod3,coeg, x xyadd,xinod1,xinod2,yinod1,yinod2,icheck) * if(icheck.eq.1) goto 50 * call trich(lhp,npoin,itot,inod1,inod2,inod3,coeg, x xyadd,xinod1,xinod2,yinod1,yinod2,icheck) * if(icheck.eq.1) goto 50 * 100 continue goto 110 50 continue write(6,*) ' Error: generation breakdown, element no.',nt ierr=9000 return 110 continue * ********************************************************* * new point found, fill in information ********************************************************* * if(inod3.gt.npoin) goto 150 * * old point, check all possibilities * call oldpo(h1,lpoin,ikno,lfapo,lface,itnode,nt, x nface,nfapo,npoin,nheap,rface,lheap,coeg) * goto 200 150 continue * ********************************************************* * fill in info about new node ********************************************************* * xp=xyadd(1,inod3-npoin) yp=xyadd(2,inod3-npoin) npoin=npoin+1 coeg(1,npoin)=xp coeg(2,npoin)=yp ibou(npoin)=-1 ikno(3)=npoin * call addq(lquad,nquad,xyquad,npoin,coeg) nodnr=inod1 call fiface(h1,nodnr,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) nodnr=inod2 call fiface(h1,nodnr,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) 200 continue nt=nt+1 * write(6,*) 'nt=',nt itnode(1,nt)=ikno(1) itnode(2,nt)=ikno(2) itnode(3,nt)=ikno(3) * if(nt.eq.23) then * write(*,*) 'nt=',nt * call dwmatg(lface,3,nface,0,'i ','lface') * stop * endif * if(nt.eq.mele) then write(6,*) ' Error: too many elements:',nt ierr=12 return endif * goto 10 9000 continue write(6,*) ' Number of generated triangles:',nt write(6,*) ' Number of generated nodes:',npoin return end * * * subroutine sort(lhp,rhp,imax) implicit double precision (a-h,o-z) dimension lhp(*),rhp(*) i=1 10 continue if(i.ge.imax) goto 9000 tal=rhp(i) lp=lhp(i) k=i+1 30 continue do 20 j=k,imax if(rhp(j).ge.tal) goto 20 lhp(i)=lhp(j) lhp(j)=lp rhp(i)=rhp(j) rhp(j)=tal goto 10 20 continue i=i+1 goto 10 9000 continue return end * * * subroutine fiface(h1,nodnr,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) * ********************************************************* * fill in info about new face ********************************************************* * implicit double precision (a-h,o-z) dimension rface(*),lheap(*),lfapo(3,*),lpoin(*),lface(3,*) dimension ikno(*),coeg(2,*),inod(2) * iface=nface+1 nface=nface+1 10 continue * ********************************************************* * nodnr is not on the face, determine nodes on face ********************************************************* * do 20 i=1,3 if(ikno(i).eq.nodnr) goto 30 20 continue 30 continue if(i.eq.1) inod(1)=ikno(3) if(i.eq.1) inod(2)=ikno(2) if(i.eq.2) inod(1)=ikno(1) if(i.eq.2) inod(2)=ikno(3) * lface(1,iface)=inod(1) lface(2,iface)=inod(2) lface(3,iface)=1 dx=coeg(1,inod(1))-coeg(1,inod(2)) dy=coeg(2,inod(1))-coeg(2,inod(2)) rface(iface)=dsqrt(dx*dx+dy*dy) call addh(lheap,rface,nheap,iface) * do 40 j=1,2 istart=lpoin(inod(j)) if(istart.eq.0) goto 50 60 continue itr=lfapo(3,istart) if(itr.lt.0) istart=-itr if(itr.lt.0) goto 60 if(itr.eq.2) lfapo(3,istart)=-(nfapo+1) if(itr.eq.2) istart=nfapo+1 if(itr.eq.2) nfapo=nfapo+1 if(itr.eq.2) itr=0 itr=itr+1 lfapo(itr,istart)=iface lfapo(3,istart)=itr goto 40 50 continue istart=nfapo+1 nfapo=nfapo+1 itr=1 lpoin(inod(j))=istart lfapo(itr,istart)=iface lfapo(3,istart)=itr 40 continue return end * * * subroutine nodcom(ikno,iface,lface,nodnr,nodco) implicit double precision (a-h,o-z) * ********************************************************* * number of nodes in common ********************************************************* * dimension ikno(*),lface(3,*) i1=0 do i=1,3 inod=ikno(i) do j=1,2 inod2=lface(j,iface) if(inod.eq.inod2) i1=i1+1 enddo enddo nodco=i1 nodnr=0 if (nodco.ne.2) goto 9000 nod3=ikno(3) nodnr=lface(1,iface) if(nodnr.eq.nod3) nodnr=lface(2,iface) 9000 continue return end subroutine oldpo(h1,lpoin,ikno,lfapo,lface,itnode,nt, x nface,nfapo,npoin,nheap,rface,lheap,coeg) implicit double precision (a-h,o-z) * dimension lface(3,*),itnode(3,*),lpoin(*),lfapo(3,*) dimension ikno(*),lheap(*),coeg(2,*),rface(*) * ********************************************************* * first: check if any old face is in new triangle ********************************************************* * iretu=0 inod1=ikno(1) inod2=ikno(2) inod3=ikno(3) istart=lpoin(inod3) 10 continue itr=lfapo(3,istart) if (itr.lt.0) itr=2 * do 30 ik=1,itr iface=lfapo(ik,istart) if(lface(3,iface).eq.0) goto 30 call nodcom(ikno,iface,lface,nodnr,nodco) * if (nodco.eq.2) iretu=iretu+1 if (nodco.eq.2.and.iretu.eq.1) nodrem=nodnr if (nodco.eq.2) goto 40 goto 30 40 continue * lfa = lface(3,iface) lface(3,iface) = lfa - 1 * * check if this boundary is to be triangulated again * if(lface(3,iface).gt.0) then do i=1,3 j=i+1-(i+1)/4*3 * if(lface(1,iface).eq.ikno(i)) then if(lface(2,iface).eq.ikno(j)) then lfa=lface(1,iface) lface(1,iface)=lface(2,iface) lface(2,iface)=lfa endif goto 30 endif enddo endif 30 continue * if (lfapo(3,istart).ge.0) go to 50 istart=-lfapo(3,istart) go to 10 50 continue * ********************************************************* * fill in info for new face (if only one old) ********************************************************* * if(iretu.eq.1) call fiface(h1,nodrem,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) if(iretu.gt.0) goto 9000 * ********************************************************* * two new faces ********************************************************* * nodnr=inod1 call fiface(h1,nodnr,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) nodnr=inod2 call fiface(h1,nodnr,lface,nface,lpoin,lfapo, x nfapo,nheap,rface,lheap,ikno,coeg,nt) 9000 continue return end subroutine poinqq(lquad,xyquad,xp,yp,radi, x iqp,npq,coeg) implicit double precision (a-h,o-z) dimension lquad(7,*),iqp(*),coeg(2,*),xyquad(3,*) * rad=radi xmax=xp+rad xmin=xp-rad ymax=yp+rad ymin=yp-rad npq=0 iquad=1 itest=1 go to 20 10 continue itest=lquad(5,iq)+1 20 continue if (itest.le.4) then iq=lquad(itest,iquad) * ********************************************************* * check limits for iq ********************************************************* * xma2=xyquad(1,iq)+xyquad(3,iq)/2.d0 xmi2=xma2-xyquad(3,iq) yma2=xyquad(2,iq)+xyquad(3,iq)/2.d0 ymi2=yma2-xyquad(3,iq) * xt1=xma2 yt1=yma2 xt2=xma2 yt2=ymi2 xt3=xmi2 yt3=yma2 xt4=xmi2 yt4=ymi2 * xma1=dmax1(xt1,xt2,xt3,xt4) xmi1=dmin1(xt1,xt2,xt3,xt4) yma1=dmax1(yt1,yt2,yt3,yt4) ymi1=dmin1(yt1,yt2,yt3,yt4) * if(ymax.ge.yma1.and.ymin.le.ymi1. $ and.xmax.ge.xma1.and.xmin.le.xmi1) go to 30 if(yma1.ge.ymax.and.ymi1.le.ymin. $ and.xma1.ge.xmax.and.xmi1.le.xmin) go to 70 if (ymax.lt.ymi1.or.ymin.gt.yma1. $ or.xmax.lt.xmi1.or.xmin.gt.xma1) go to 10 go to 70 * ********************************************************* * itest=5, up one level ********************************************************* * else if (iquad.eq.1) go to 9000 itest=lquad(5,iquad)+1 iquad=lquad(6,iquad) go to 20 endif 30 continue * ********************************************************* * iq within limits, take all points ********************************************************* * itest=1 iqold=iquad iquad=iq 40 continue * ********************************************************* * look down ********************************************************* * ip=lquad(7,iquad) * ********************************************************* * ip.ge.0 fill in points ********************************************************* * if (ip.ge.0) then if (ip.eq.0) go to 60 do ijk=1,ip iqp(npq+ijk)=lquad(ijk,iquad) enddo npq=npq+ip 60 continue * ********************************************************* * next quad ********************************************************* * itest=lquad(5,iquad)+1 iquad=lquad(6,iquad) if (iquad.eq.iqold) go to 10 if (itest.eq.5) go to 60 go to 40 else * ********************************************************* * ip.lt.0 go down ********************************************************* * iquad=lquad(itest,iquad) itest=1 go to 40 endif 70 continue * ********************************************************* * part of iq within limits, look down ********************************************************* * ip=lquad(7,iq) if (ip.lt.0) then * ********************************************************* * has sons, down one level ********************************************************* * iquad=iq itest=1 go to 20 else * ********************************************************* * no sons, check all points ********************************************************* * if (ip.eq.0) go to 10 do 80 ipo=1,ip ipt=lquad(ipo,iq) xpt=coeg(1,ipt) ypt=coeg(2,ipt) dxt=xpt-xp dyt=ypt-yp ra=dsqrt(dxt*dxt+dyt*dyt) if (radi.lt.ra) go to 80 npq=npq+1 iqp(npq)=lquad(ipo,iq) 80 continue go to 10 endif 9000 continue return end subroutine fcross(xm,ym,iv3,npoin,coeg,iface,lface, x icheck,xyadd) implicit double precision (a-h,o-z) dimension coeg(2,*),lface(3,*),xyadd(2,*) * ********************************************************* * do the segments cross ? * icheck=0 no * icheck=1 yes ********************************************************* * eps=1.d-9 icheck=0 * ********************************************************* * find point of intersection ********************************************************* * if((iv3-npoin).gt.0) then xiv3=xyadd(1,iv3-npoin) yiv3=xyadd(2,iv3-npoin) else xiv3=coeg(1,iv3) yiv3=coeg(2,iv3) endif * iv1=lface(1,iface) iv2=lface(2,iface) dx1=xiv3-xm dy1=yiv3-ym xiv1=coeg(1,iv1) yiv1=coeg(2,iv1) xiv2=coeg(1,iv2) yiv2=coeg(2,iv2) dx2=xiv2-xiv1 dy2=yiv2-yiv1 det=dx2*dy1-dy2*dx1 if(abs(det).le.eps) goto 9000 scut=(dy2*(xm-xiv1)+dx2*(yiv1-ym))/det tcut=(dy1*(xm-xiv1)+dx1*(yiv1-ym))/det * if(scut.gt.0.d0.and.scut.lt.1.d0.and. x tcut.gt.0.d0.and.tcut.lt.1.d0) icheck=1 * 9000 continue return end * * * subroutine relist(lhp,rhp,xp,yp,coeg,h1,imax) implicit double precision (a-h,o-z) dimension lhp(*),coeg(2,*),rhp(*) * ********************************************************* * check if new point is to be on top of list ********************************************************* * if(imax.eq.5) goto 9000 iglob=lhp(imax-4) i1=imax-5 h1=h1 do 35 j=1,i1 i2=lhp(j) xi2=coeg(1,i2) yi2=coeg(2,i2) dx1=xp-xi2 dy1=yp-yi2 r1=dsqrt(dx1*dx1+dy1*dy1) * if(r1.ge.0.6d0*h1) then do 40 i=j,i1 is=i1+j-i lhp(is+1)=lhp(is) rhp(is+1)=rhp(is) 40 continue lhp(j)=iglob goto 9000 endif 35 continue 9000 continue return end * * * subroutine facech(lhp,npoin,imax,lface,lpoin,lfapo,inod3,coeg, x xyadd,xinod1,xinod2,yinod1,yinod2,icheck) implicit double precision(a-h,o-z) dimension lhp(*),lfapo(3,*),lpoin(*),coeg(2,*),lface(3,*) dimension xyadd(2,*) * ********************************************************* * check that no face (adjacent to * the points) is crossed ********************************************************* * icheck=0 do 80 ij=1,imax if(lhp(ij).gt.npoin) goto 80 istart=lpoin(lhp(ij)) 90 continue i2=0 itr=lfapo(3,istart) if(itr.lt.0) itr=2 do 100 ik=1,itr iface=lfapo(ik,istart) * ********************************************************* * node on face ? ********************************************************* * if(lface(3,iface).eq.0) goto 100 * if(lface(1,iface).eq.inod3.or.lface(2,iface). x eq.inod3) goto 100 * call fcross(xinod1,yinod1,inod3,npoin,coeg,iface,lface, x icheck, xyadd) if(icheck.eq.1) goto 50 call fcross(xinod2,yinod2,inod3,npoin,coeg,iface,lface, x icheck, xyadd) if(icheck.eq.1) goto 50 100 continue if(lfapo(3,istart).lt.0) i2=1 if(i2.eq.0) goto 80 istart=-lfapo(3,istart) goto 90 80 continue goto 110 50 continue 110 continue return end * * * subroutine ficlos(lquad,x,y,xyquad,inr,nrnr) implicit double precision (a-h,o-z) dimension xyquad(3,*),lquad(7,*),inr(*) * nrnr=0 iquad=1 call findp(lquad,x,y,iquad,xyquad) * iquad=iquad if(lquad(7,iquad).eq.0) iquad=lquad(6,iquad) call pinq(inr,lquad,iquad,nrnr) * return end * * * subroutine pinq(iqp,lquad,iquadt,ippp) implicit double precision (a-h,o-z) dimension lquad(7,*),iqp(*) * npq=lquad(7,iquadt) if (npq.gt.0) then do 5 ipq=1,npq ippp=ippp+1 iqp(ippp)=lquad(ipq,iquadt) 5 continue return endif * ********************************************************* * take all points in iquadt ********************************************************* * itest=1 iqold=iquadt iquad=lquad(itest,iquadt) 40 continue * ********************************************************* * look down ********************************************************* * ip=lquad(7,iquad) * ********************************************************* * ip.ge.0 fill in points ********************************************************* * if (ip.ge.0) then if (ip.eq.0) go to 60 do ijk=1,ip iqp(ippp+ijk)=lquad(ijk,iquad) enddo ippp=ippp+ip 60 continue * ********************************************************* * back one level ********************************************************* * itest=lquad(5,iquad)+1 iquad=lquad(6,iquad) if (iquad.eq.iqold) return if (itest.eq.5) goto 60 go to 40 else * ********************************************************* * ip.lt.0 go down ********************************************************* * iquad=lquad(itest,iquad) itest=1 go to 40 endif return end * * * subroutine scdef(pu,pv,vx,vy,sc) implicit double precision (a-h,o-z) dimension pu(*),pv(*) * sc=0.d0 do 10 i=1,3 sc=sc+pu(i)*pu(i)*vx*vx sc=sc+2.d0*pu(i)*pv(i)*vx*vy sc=sc+pv(i)*pv(i)*vy*vy 10 continue * sc=dsqrt(sc) * return end * * * subroutine trich(lhp,npoin,imax,inod1,inod2,inod3,coeg, x xyadd,xinod1,xinod2,yinod1,yinod2,icheck) implicit double precision(a-h,o-z) dimension lhp(*),coeg(2,*),xyadd(2,*),c(3) * ********************************************************* * check that the triangle does not cointain a point ********************************************************* * icheck=0 if((inod3-npoin).gt.0) then xiv3=xyadd(1,inod3-npoin) yiv3=xyadd(2,inod3-npoin) else xiv3=coeg(1,inod3) yiv3=coeg(2,inod3) endif * do 80 ij=1,imax iact=lhp(ij) if(iact.gt.npoin) goto 80 * if(iact.eq.inod1.or.iact.eq.inod2.or.iact.eq.inod3) goto 80 * ********************************************************* * check barycentric coordinates ********************************************************* * x=coeg(1,iact) y=coeg(2,iact) * x2=xinod2-xinod1 y2=yinod2-yinod1 x3=xiv3-xinod1 y3=yiv3-yinod1 xr=x-xinod1 yr=y-yinod1 det=x2*y3-x3*y2 c(2)=(xr*y3-x3*yr)/det c(3)=(x2*yr-xr*y2)/det c(1)=1.0d0-c(2)-c(3) * i1=0 do 70 i=1,3 if(c(i).ge.0.d0.and.c(i).le.1.d0) i1=i1+1 70 continue if(i1.eq.3) then icheck=1 return endif * 80 continue return end * * * subroutine arheap( i ifarea,lface,coeg, o lheap,rface,nheap) implicit double precision(a-h,o-z) dimension lface(3,*),coeg(2,*) dimension ifarea(*),lheap(*),rface(*) * ********************************************************* * adds faces on current area to the heap * * ifarea(*) faces on areas * (1): number of faces on current area* * (2-): face numbers * lheap(nheap) heaplist for segments * rface(nheap) segment sizes ********************************************************* * i1=0 nface=ifarea(1) nmax=nface+1 * do 20 ij=2,nmax i=ifarea(ij) * ip1=lface(1,i) ip2=lface(2,i) dx=coeg(1,ip1)-coeg(1,ip2) dy=coeg(2,ip1)-coeg(2,ip2) dl=dx*dx+dy*dy i1=i1+1 rface(i1)=dsqrt(dl) if(i1.eq.1) then lheap(1)=i nheap=1 else call addh(lheap,rface,nheap,i) endif 20 continue return end * * * subroutine lfapc( i lface,nface,npoin, o lpoin,lfapo,nmax) implicit double precision(a-h,o-z) dimension lface(3,*),lpoin(*),lfapo(3,*) * ********************************************************* * lpoin(nno) connection node nr - position in lfapo * lfapo(3,*) linked list of node-face couplings * 1-2: > 0 - a face surrounding a node * 3: number of faces stored at * this location (1 or 2) ********************************************************* * nmax=npoin+1 * do 300 ij=1,npoin lpoin(ij)=ij do 200 ipp=1,nface do 100 j=1,2 is=ij if(lface(j,ipp).eq.is) then 50 continue itr=lfapo(3,is) if(itr.lt.0) then is=-itr goto 50 endif if(itr.lt.2) then lfapo(itr+1,is)=ipp lfapo(3,is)=lfapo(3,is)+1 endif if(itr.eq.2) then lfapo(3,is)=-nmax lfapo(1,nmax)=ipp lfapo(3,nmax)=1 nmax=nmax+1 endif endif 100 continue 200 continue 300 continue return end subroutine straight(ipoin,ip1,ip2,coeg,x,y,h1) * ********************************************************* * coeg=coordinates for points * * x,y input old location, output new location ********************************************************* * implicit double precision (a-h,o-z) dimension coeg(2,*) * dx=coeg(1,ip2)-coeg(1,ip1) dy=coeg(2,ip2)-coeg(2,ip1) dltot=dsqrt(dx*dx+dy*dy) * dx=coeg(1,ip2)-coeg(1,ipoin) dy=coeg(2,ip2)-coeg(2,ipoin) * dl=dsqrt(dx*dx+dy*dy) vx=dx/dl vy=dy/dl * xp=x yp=y * xnew=h1*vx+xp ynew=h1*vy+yp * xt=coeg(1,ip2)-xnew yt=coeg(2,ip2)-ynew dleft=dsqrt(xt*xt+yt*yt) dlx=xnew-coeg(1,ip1) dly=ynew-coeg(2,ip1) dll=dsqrt(dlx*dlx+dly*dly) if(dleft.lt.0.5d0*h1.or.dll.ge.dltot) then x=coeg(1,ip2) y=coeg(2,ip2) ipoin=ip2 return endif x=xnew y=ynew return end subroutine spline(ipoin,ipstep,ip1,ip2,coeg,told, $ x,y,vn,h1,nextra,cseg,nrcseg,inamcs) * ********************************************************* * coeg=coordinates for points * * x,y input old location, output new location * * ipstep shows how far we have gone along the spline ********************************************************* * implicit none integer nrcseg,inamcs,ipoin,ipstep,ip1,ip2,ipstepin,nextra,i double precision coeg(2,*),cseg(nrcseg,4,*) double precision p1(2),p2(2),v1(2),v2(2),p(2),vn(2) double precision dx,dy,dltot,h1,hleft,dl,dl1,dlen,dleft double precision x,y,xnew,ynew,vx,vy,t,hact,hin,told * dx=coeg(1,ip2)-coeg(1,ip1) dy=coeg(2,ip2)-coeg(2,ip1) dltot=dsqrt(dx*dx+dy*dy) * ********************************************************* * normal vectors and coordinates ********************************************************* * ipstepin=ipstep hin=h1 hleft=h1 10 continue if(ipstep.eq.1) then p1(1)=coeg(1,ip1) p1(2)=coeg(2,ip1) v1(1)=cseg(inamcs,1,1) v1(2)=cseg(inamcs,2,1) else p1(1)=cseg(inamcs,1,ipstep) p1(2)=cseg(inamcs,2,ipstep) v1(1)=cseg(inamcs,3,ipstep) v1(2)=cseg(inamcs,4,ipstep) endif * if(ipstep.eq.nextra+1) then p2(1)=coeg(1,ip2) p2(2)=coeg(2,ip2) v2(1)=cseg(inamcs,3,1) v2(2)=cseg(inamcs,4,1) else p2(1)=cseg(inamcs,1,ipstep+1) p2(2)=cseg(inamcs,2,ipstep+1) v2(1)=cseg(inamcs,3,ipstep+1) v2(2)=cseg(inamcs,4,ipstep+1) endif * ********************************************************* * try out the location of the new point ********************************************************* * 20 continue dx=p2(1)-coeg(1,ipoin) dy=p2(2)-coeg(2,ipoin) dl=dsqrt(dx*dx+dy*dy) vx=dx/dl vy=dy/dl xnew=h1*vx+x ynew=h1*vy+y dl1=dsqrt((xnew-x)*(xnew-x)+(ynew-y)*(ynew-y)) * ********************************************************* * do we have to pass by a control point? ********************************************************* * if(dl1.gt.dl) then if(ipstep.eq.nextra+1) goto 100 hleft=dmax1(hleft-dl,0.d0) ipstep=ipstep+1 told=0.d0 goto 10 endif * ********************************************************* * find actual position using hermitian interpolation ********************************************************* * dx=p2(1)-p1(1) dy=p2(2)-p1(2) dlen=dsqrt(dx*dx+dy*dy) t=told+hleft/dlen write(6,*) 'hleft=',hleft,' dlen=',dlen write(6,*) 'hleft/dlen=',hleft/dlen write(6,*) 't=',t,' told=',told call findspl(t,p1,p2,v1,v2,p,vn) xnew=p(1) ynew=p(2) write(6,*) 'p1=',(p1(i),i=1,2) write(6,*) 'p2=',(p2(i),i=1,2) write(6,*) 'v1=',(v1(i),i=1,2) write(6,*) 'v2=',(v2(i),i=1,2) write(6,*) 'xnew2=',xnew,' ynew2=',ynew write(6,*) 'xold=',coeg(1,ipoin),' yold=',coeg(2,ipoin) hact=dsqrt((xnew-x)*(xnew-x)+(ynew-y)*(ynew-y)) * ********************************************************* * is the actual h close enough to the desired h? ********************************************************* * write(6,*) 'hact=',hact,' h1=',h1 write(6,*) '******************************' if(hact.gt.1.2d0*hin.or.hact.lt.0.8d0*hin) then if(hact.gt.1.2d0*h1) then h1=0.9d0*h1 else h1=1.1d0*h1 endif ipstep=ipstepin hleft=h1 if(hleft.lt.1.d-12) stop goto 10 endif write(6,*) 'PASSED' * ********************************************************* * right meshsize, create new point ********************************************************* * told=t write(6,*) 'told=',told dx=coeg(1,ip2)-xnew dy=coeg(2,ip2)-ynew dleft=dsqrt(dx*dx+dy*dy) dx=xnew-coeg(1,ip1) dy=ynew-coeg(2,ip1) dl=dsqrt(dx*dx+dy*dy) if(dleft.lt.0.5d0*h1.or.dl.ge.dltot) goto 100 x=xnew y=ynew goto 9000 100 continue x=coeg(1,ip2) y=coeg(2,ip2) ipoin=ip2 9000 continue return end subroutine ichecc(nrcseg,icseg,isname,itrue,i) implicit none integer itrue,nrcseg,isname,icseg(4,*),i * itrue=0 do i=1,nrcseg if(icseg(1,i).eq.isname) then itrue=1 return endif enddo return end subroutine shuffle(cseg,nrcseg,inamcs,nextra) implicit none integer nrcseg,inamcs,nextra,i,j,k,l double precision cseg(nrcseg,4,*),work * do i=1,2 work=cseg(inamcs,i,1) cseg(inamcs,i,1)=cseg(inamcs,i+2,1) cseg(inamcs,i+2,1)=work enddo l=(nextra+1)/2 do j=2,l k=nextra-j+3 do i=1,4 work=cseg(inamcs,i,j) cseg(inamcs,i,j)=cseg(inamcs,i,k) cseg(inamcs,i,k)=work enddo enddo return end subroutine findspl( i t,p1,p2,v1,v2, o p,v) implicit none integer i double precision t,p2(*),p1(*),v2(*),v1(*),p(*) double precision phi(2),psi(2),d2(2),d1(2),rot,v(2) * * p1,p2 endpoints * v1,v2 normal vectors at the endpoints * t t is in (0,1) in the interval (p1,p2) * p is the point above t * v is the normal vector at p * phi(1)=1.d0+(2.d0*t-3.d0)*t*t phi(2)=(3.d0-2.d0*t)*t*t psi(1)=(1.d0+(t-2.d0)*t)*t psi(2)=(t-1.d0)*t*t * rot = v1(1)*(p2(2)-p1(2))-v1(2)*(p2(1)-p1(1)) d1(1)=-rot*v1(2) d1(2)=rot*v1(1) rot = v2(1)*(p2(2)-p1(2))-v2(2)*(p2(1)-p1(1)) d2(1)=-rot*v2(2) d2(2)=rot*v2(1) * rot=dsqrt(d1(1)*d1(1)+d1(2)*d1(2)) d1(1)=d1(1)/rot d1(2)=d1(2)/rot rot=dsqrt(d2(1)*d2(1)+d2(2)*d2(2)) d2(1)=d2(1)/rot d2(2)=d2(2)/rot * do i=1,2 p(i)=phi(1)*p1(i)+phi(2)*p2(i)+psi(1)*d1(i)+psi(2)*d2(i) v(i)=(1.d0-t)*v1(i)+t*v2(i) enddo return end