program skew_mesh parameter( mx=500 ) common h0 real x_grid(mx*mx), y_grid(mx*mx), x_temp(mx), y_temp(mx) integer bb_x, bb_y character*9 moveto, lineto open(unit=1,file='mesh.inp') open(unit=4,file='grid.dat') c open(unit=7,file='mesh.ps') open(unit=7,file='mesh.eps') c *** Read input parameters. read(1,5 ) nx, ny, h0, sf_x, sf_y, tol 5 format( i10,/,i10,/,f10.5,/,f10.5,/,f10.5,/,e10.1 ) if( nx .gt. mx .or. ny .gt. mx ) then print *, 'nx or ny is too big, enlarge parameter mx' stop end if c *** Determine the grid point distirbution on the perimeter. call outline( nx, ny, sf_x, sf_y, x_grid, y_grid ) c *** Generate an initial mesh on the interior. call initial_grid( nx, ny, x_grid, y_grid ) c *** Flip the y coordinates h = height(0.0) do i=1, Nx do j=1, Ny ind = Nx*(j-1)+i y_temp(j) = y_grid(ind) x_temp(j) = x_grid(ind) end do do j=1, Ny ind = Nx*(j-1)+i y_grid(ind) = h - y_temp(Ny-j+1) x_grid(ind) = x_temp(Ny-j+1) end do end do c *** Write the mesh to file mesh.dat write(4,10) nx, ny 10 format(2i10) do j=1, ny i_skip = (j-1)*nx do i=1, nx write(4,20) x_grid(i_skip+i), y_grid(i_skip+i) 20 format(1p2e15.7) end do end do c *** Write a postscript version of the mesh to file mesh.ps moveto = ' moveto' lineto = ' lineto' scale = 1.0 x_0 = x_grid(1) y_0 = y_grid(1) x_1 = x_grid(nx*ny) y_1 = y_grid(nx*ny) x_L = x_1 - x_0 y_L = y_1 - y_0 x_fit = (11.0-1.0)*72/x_L y_fit = ( 8.5-1.0)*72/y_L scale = scale*min( x_fit, y_fit ) bb_x = int(x_L*scale) bb_y = int(y_L*scale) translate_x = ( 8.5*72 - bb_y )/2 translate_y = ( 11.0*72 - bb_x )/2 c write(7,50) '%!PS-Adobe-2.0', c & '%%BoundingBox: ',0, 0, bb_x, bb_y, c & translate_x, translate_y, ' translate', c & '90 rotate', c & '1 -1 scale', c & '0.2 setlinewidth' c50 format(a14,/,a15,4i5,//,2f8.2,a11,/a11,/,a12,//,a16,/) margin = 10 write(7,50) '%!PS-Adobe-2.0 EPSF-2.0', & '%%BoundingBox: ', & -margin, -margin, bb_x+margin, bb_y+margin, & '0.2 setlinewidth' 50 format(a23,/,a15,4i5,//,a16,/) do j=1, ny ind_0 = nx*(j-1) write(7,70) scale*(x_grid(ind_0+1)-x_0), & scale*(y_grid(ind_0+1)-y_0), moveto 70 format(2f10.4,a9) do i=2, nx ind = ind_0 + i write(7,70) scale*(x_grid(ind)-x_0), & scale*(y_grid(ind)-y_0), lineto end do write(7,80) 'stroke' 80 format(a6) end do do i=1, nx write(7,70) scale*(x_grid(i)-x_0), & scale*(y_grid(i)-y_0), moveto do j=2, ny ind = nx*(j-1) + i write(7,70) scale*(x_grid(ind)-x_0), & scale*(y_grid(ind)-y_0), lineto end do write(7,80) 'stroke' end do write(7,90) 'showpage' 90 format(a8) stop end c------------------------------------------------------------------------------ subroutine initial_grid( nx, ny, x_grid, y_grid ) dimension x_grid(nx,ny), y_grid(nx,ny) c *** Average the mesh from the boundaries to the interior yL_1 = 1.0/( y_grid(1,ny) - y_grid(1,1) ) do j=2, ny-1 w1 = ( y_grid(1,j) - y_grid(1,1) )*yL_1 w2 = 1.0 - w1 do i=2, nx-1 x_grid(i,j) = w2*x_grid(i,1) + w1*x_grid(i,ny) y_grid(i,j) = w2*y_grid(i,1) + w1*y_grid(i,ny) end do end do return end c------------------------------------------------------------------------------ subroutine outline( nx, ny, sf_x, sf_y, x_grid, y_grid ) parameter( mx=500 ) dimension x_grid(nx,ny), y_grid(nx,ny) dx = 1.0/float(Nx-1) dy = height(0.0)/float(Ny-1) do i=(Nx+1)/2, Nx isymm = Nx+1-i x = htgrid(0.5,1.0,i-(Nx+1)/2+1,(Nx+1)/2,sf_x) x_grid(i ,1) = x x_grid(isymm,1) = 1.0 - x y_grid(i ,1) = height(0.0)-height(x) y_grid(isymm,1) = height(0.0)-height(x) x_grid(i ,Ny) = x x_grid(isymm,Ny) = 1.0 - x y_grid(i ,Ny) = height(0.0) y_grid(isymm,Ny) = height(0.0) end do do j=2, Ny-1 y = htgrid(y_grid(1,j),y_grid(1,Ny),j,Ny,sf_y) x_grid(1,j) = 0 y_grid(1,j) = y x_grid(Nx,j) = 1.0 y_grid(Nx,j) = y end do return end c------------------------------------------------------------------------------ function height(x) common h0 pi = acos(-1.0) height = h0*( 2.0 - sin( pi*x )**4 ) return end c------------------------------------------------------------------------------ function htgrid(x0,xL,i,N,a) if( a .eq. 0.0 ) then htgrid = float(i-1)/float(N-1)*( xL - x0 ) + x0 else atanha=0.5*log(( 1.0+a)/(1.0-a)) xi = 1-float(i-1)/float(N-1) htgrid = ( 1 - tanh(xi*atanha)/a ) * ( xL - x0 ) + x0 end if return end c------------------------------------------------------------------------------ function htfact(y1,dy1,ym,m) tol=1.e-8 itrmax=20 N = m if( dy1 .gt. (ym-y1)/(m-1) ) then print *, 'ERROR in htfact: dy1 .gt. (ym-y1)/(m-1)' stop end if xi_2 = 1 - 1.0/(N-1) xi_m = 1 - float(m-1)/(N-1) ratio = dy1/(ym-y1) a = 0.95 do i=1,itrmax atanha = 0.5*log( (1.0+a)/(1.0-a) ) a_inv = 1.0 / a f2 = 1 - tanh(xi_2*atanha)*a_inv fm = 1 - tanh(xi_m*atanha)*a_inv df2 = ( 1-f2 - xi_2/( (1-a*a)*cosh(xi_2*atanha)**2 ) )*a_inv dfm = ( 1-fm - xi_m/( (1-a*a)*cosh(xi_m*atanha)**2 ) )*a_inv delta = ( ratio*fm - f2 ) / ( df2 - f2/fm*dfm ) a = a + delta if( a .gt. 1 ) a = a - 0.5*delta if( abs(delta/a) .lt. tol ) goto 5 end do print *, & 'WARNING from htfact: iterations did not converge, delta = ', & delta 5 continue htfact = a return end