パリティ Vol.06 No.03 (1991) 62. のフォートラン入門より、平行平板コンデンサの作る二次元電場の電位を、ラプラス方程式のヤコビによる解法で解きます。
forallを用い mask 配列を使ってコンデンサの電極は計算を飛ばすことにします。forall は並列代入操作で、ただの do loop の時と配列更新のタイミングが違い、値は全ループが終わるまで更新されません。
さぼって収束性は確かめずに、十分と思われる回数反復することにします。計算の核になる部分はほんとに短く書けます。
ラプラス方程式は、熱伝導やゴムを伸ばした時の形状も記述しているのでコンデンサの電極の代わりに自由に定値を与えてやれば、ラバーフェチも満足の計算ができます。
グラフィックは、OpenGLでイタレーションごとに描きたかったのですが、めんどくさかったので収束後の電位だけを、いい加減にプロットすることにしましたw
ソース・プログラム
program jacobi use plotter implicit none integer, parameter :: nsize = 50, nsz = nsize + 1 real :: v(-nsz:nsz, -nsz:nsz) = 0.0 logical :: mask(-nsz:nsz, -nsz:nsz) = .true. integer :: i, j, iter mask(-10, -5:5) = .false. ! static voltage mask( 10, -5:5) = .false. v(-10, -5:5) = 10.0 v( 10, -5:5) = -10.0 ! do iter = 1, 200 forall(i = -nsize:nsize, j = -nsize:nsize, mask(i, j)) & v(i, j) = 0.25 * ( v(i - 1, j) + v(i + 1, j) + v(i, j - 1) + v(i, j + 1) ) end do ! call graph() stop contains subroutine graph() integer :: ix, iy real :: x, y call gr_on('Laplace equation', 600, 600) call gr_pencol(0) ! x-direction do iy = -nsize, nsize, 2 x = 0.5 + 0.005 * -nsize + 0.005 * iy y = 0.5 + 0.003 * iy call move( x, y ) do ix = -nsize, nsize, 2 x = 0.005 * ix + 0.005 * iy + 0.5 y = 0.5 + 0.003 * iy + v(ix, iy) * 0.02 call draw( x, y ) end do end do ! y-direction do ix = -nsize, nsize, 2 x = 0.005 * ix + 0.005 * -nsize + 0.5 y = 0.5 + 0.003 * -nsize call move( x, y ) do iy = -nsize, nsize, 2 x = 0.005 * ix + 0.005 * iy + 0.5 y = 0.5 + 0.003 * iy + v(ix, iy) * 0.02 call draw( x, y ) end do end do call gr_off() return end subroutine graph end program jacobi
グラフィックス・ルーチン
module uho_win use ifwina use ifwinty use ifmt, only : RTL_CRITICAL_SECTION implicit none type :: t_wnd integer (HANDLE) :: hWnd integer (HANDLE) :: hDC integer (LPINT) :: hThread integer :: id integer (HANDLE) :: hPen character (LEN = 80) :: title = 'Fortran Plot' integer :: nsize_x = 320, nsize_y = 240 real :: xmin = 0.0, xmax = 1.0, ymin = 0.0, ymax = 1.0 end type type (t_wnd) :: wnd type (RTL_CRITICAL_SECTION) :: lpCriticalSection contains !-------------------------------------------------------------------------------- integer(4) function WinMain( hInstance, hPrevInstance, lpszCmdLine, nCmdShow, wnd) implicit none integer (HANDLE), intent(in) :: hInstance, hPrevInstance integer (LPSTR) , intent(in) :: lpszCmdLine integer (SINT) , intent(in) :: nCmdShow type (t_wnd) , intent(in out) :: wnd ! Variables logical, save :: qfirst = .true. type (T_WNDCLASS) :: wc type (T_MSG) :: mesg integer (HANDLE) :: hWndMain integer (BOOL) :: iretb character (LEN = 256) :: ClassName = 'Fortran'//char(0) integer :: noffset_x, noffset_y !side line = 6, title bar = 25 ! Init Main window noffset_x = 2 * GetSystemMetrics(SM_CXFIXEDFRAME) noffset_y = 2 * GetSystemMetrics(SM_CYFIXEDFRAME) + GetSystemMetrics(SM_CYCAPTION) ! WinMain = -1 ! Error code if (qfirst) then wc%lpszClassName = loc(ClassName) wc%lpfnWndProc = loc(MainWndProc) ! call BACK procedure name wc%style = ior(CS_VREDRAW , CS_HREDRAW) wc%hInstance = hInstance wc%hIcon = NULL wc%hCursor = LoadCursor( NULL, IDC_ARROW ) wc%hbrBackground = ( COLOR_WINDOW + 1 ) if ( RegisterClass(wc) == 0 ) return qfirst = .FALSE. end if !Init instance WinMain = -2 ! Error code hWndMain = CreateWindowEx( 0, ClassName, & trim(wnd%title)//char(0), & int(ior(WS_OVERLAPPED, WS_SYSMENU)), & CW_USEDEFAULT, CW_USEDEFAULT, & wnd%nsize_x + noffset_x, & wnd%nsize_y + noffset_y, & 0, & 0, & hInstance, & loc(wnd) ) if (hWndMain == 0) return iretb = ShowWindow( hWndMain, nCmdShow ) iretb = UpdateWindow( hWndMain ) ! Message Loop do while ( GetMessage (mesg, NULL, 0, 0) ) iretb = TranslateMessage( mesg ) iretb = DispatchMessage( mesg ) end do WinMain = mesg%wParam return end function WinMain !---------------------------------------------------------- integer (LRESULT) function MainWndProc( hWnd, mesg, wParam, lParam ) !DEC$ ATTRIBUTES STDcall, DECORATE, ALIAS : 'MainWndProc' :: MainWndProc use IFMT implicit none integer (HANDLE) , intent(in) :: hWnd integer (UINT) , intent(in) :: mesg integer (fwParam), intent(in) :: wParam integer (flParam), intent(in) :: lParam ! integer (HANDLE) :: hDC, hBmp integer (BOOL) :: iretb type (T_PAINTSTRUCT) :: ps type (T_RECT) :: rc ! MainWndProc = 0 select case ( mesg ) case (WM_CREATE) hDC = GetDC(hWnd) iretb = GetClientRect(hWnd, rc) hBmp = CreateCompatibleBitmap(hDC, rc%right - rc%left, rc%bottom - rc%top) wnd%hWnd = hWnd wnd%hDC = CreateCompatibleDC(hDC) iretb = SelectObject(wnd%hDC, hBmp) iretb = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS) iretb = ReleaseDC(hWnd, hDC) iretb = DeleteObject(hBmp) case (WM_DESTROY) iretb = DeleteObject( wnd%hDC ) call PostQuitMessage( 0 ) case (WM_PAINT) call EnterCriticalSection( loc(lpCriticalSection) ) hDC = BeginPaint( wnd%hWnd, ps ) iretb = GetClientRect( wnd%hWnd, rc ) iretb = BitBlt(hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, wnd%hDC, 0, 0, SRCCOPY) iretb = endPaint( wnd%hWnd, ps ) call LeaveCriticalSection( loc(lpCriticalSection) ) case (WM_LBUTTONDOWN) write(*, *) 'LBUTTON ', trim(wnd%title), wnd%hWnd case (WM_RBUTTONDOWN) iretb = DeleteObject( wnd%hDC ) call PostQuitMessage( 0 ) case default MainWndProc = DefWindowProc( hWnd, mesg, wParam, lParam ) end select return end function MainWndProc !------------------------------------------------------------------------------- integer (HANDLE) function get_DosHndl() implicit none interface function GetConsoleWindow() ! non-existent in IFORTxx.MOD !DEC$ ATTRIBUTES DEFAULT, STDcall, DECORATE, ALIAS:'GetConsoleWindow' :: GetConsoleWindow use ifwinty integer (HANDLE) :: GetConsoleWindow end function end interface get_DosHndl = GetConsoleWindow() return end function get_DosHndl !------------------------------------------------------------------------------------- end module uho_win !================================================================================= module uhoplot use uho_win contains !------------------------------------------------------------------------------------- subroutine gr_on(text, nx, ny) use IFMT ! multithread module implicit none character (LEN = *), intent(in), optional :: text integer , intent(in), optional :: nx, ny integer (BOOL) :: iretb integer (HANDLE) :: hBmp integer (LPDWORD) :: id type (T_RECT) :: rc if ( present(text) ) wnd%title = TRIM(text) if ( present(nx) ) wnd%nsize_x = nx if ( present(ny) ) wnd%nsize_y = ny call InitializeCriticalSection( loc(lpCriticalSection)) wnd%hThread = CreateThread(0, 0, Thread_Proc, loc(wnd), CREATE_SUSPendED, id) wnd%id = id iretb = SetThreadPriority(wnd%hThread, THREAD_PRIORITY_BELOW_NORMAL) iretb = ResumeThread(wnd%hThread) call sleep(100) ! wait for Window initialization iretb = GetClientRect(wnd%hWnd, rc) hBmp = CreateCompatibleBitmap(wnd%hDC, rc%right - rc%left, rc%bottom - rc%top) iretb = SelectObject(wnd%hDC, hBmp) iretb = DeleteObject(hBmp) iretb = PatBlt(wnd%hDC, 0, 0, rc%right - rc%left, rc%bottom - rc%top, WHITENESS) wnd%hPen = CreatePen(PS_SOLID, 1, 0) return end subroutine gr_on !------------------------------------------------------------------------------------- subroutine gr_off(isec) use IFMT ! Module for multithread implicit none integer (DWORD), intent(in), optional :: isec integer (BOOL) :: iretb integer (DWORD) :: iwait if (present(isec)) THEN iwait = isec * 1000 else iwait = INFINITE end if call gr_show() iretb = DeleteObject(wnd%hPen) iretb = WaitForSingleObject(wnd%hThread, iwait) iretb = CloseHandle(wnd%hThread) iretb = PostMessage(wnd%hWnd, WM_DESTROY, NULL, NULL) wnd%hThread = NULL call DeleteCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_off !------------------------------------------------------------------------------------- integer (LONG) function Thread_Proc(lp_ThreadParameter) !DEC$ ATTRIBUTES STDcall, ALIAS:"_thread_proc" :: Thread_Proc implicit none integer (LPINT), intent(in) :: lp_ThreadParameter integer :: hInst type (t_wnd) :: wnd pointer (p_wnd, wnd) ! non-standard fortran p_wnd = lp_ThreadParameter hInst = GetModuleHandle(NULL) Thread_Proc = WinMain(hInst, NULL, NULL, SW_SHOWNORMAL, wnd) return end function Thread_Proc !------------------------------------------------------------------------------------- subroutine gr_dot(ix, iy, icol) implicit none integer, intent(in) :: ix, iy integer, intent(in) :: icol integer (BOOL):: iretb call EnterCriticalSection( loc(lpCriticalSection) ) iretb = SetPixel(wnd%hDC, ix, iy, icol) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_dot !------------------------------------------------------------------------------------- subroutine gr_move(ix, iy) implicit none integer, intent(in) :: ix, iy integer (BOOL):: iretb call EnterCriticalSection( loc(lpCriticalSection) ) iretb = MoveToEx(wnd%hDC, ix, iy, NULL) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_move !------------------------------------------------------------------------------------- subroutine gr_line(ix, iy) implicit none integer, intent(in) :: ix, iy integer (BOOL):: iretb call EnterCriticalSection( loc(lpCriticalSection) ) iretb = LineTo(wnd%hDC, ix, iy) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_line !------------------------------------------------------------------------------------- subroutine gr_show() implicit none integer (BOOL):: iretb call EnterCriticalSection( loc(lpCriticalSection) ) iretb = InvalidateRect(wnd%hWnd, NULL, FALSE) call LeaveCriticalSection( loc(lpCriticalSection)) return end subroutine gr_show !------------------------------------------------------------------------------------- subroutine gr_pencol(icol) implicit none integer, intent(in) :: icol integer :: iretb call EnterCriticalSection( loc(lpCriticalSection) ) iretb = DeleteObject(wnd%hPen) wnd%hPen = CreatePen(PS_SOLID, 1, icol) iretb = SelectObject(wnd%hDC, wnd%hPen) call LeaveCriticalSection( loc(lpCriticalSection)) return end subroutine gr_pencol !------------------------------------------------------------------------------------- subroutine gr_text(ix, iy, txt, icol, ifontsize, ifontdirection) implicit none integer, intent(in) :: ix, iy character (LEN = *), intent(in) :: txt integer, intent(in), optional :: icol integer (BOOL) :: iretb integer (HANDLE) :: hFont integer, optional, intent(in) :: ifontsize, ifontdirection integer :: kfontsize, kfontdirection if ( present(icol) ) iretb = SetTextColor(wnd%hDC, icol) if ( present(ifontsize) ) then kfontsize = ifontsize else kfontsize = 10 end if if ( present(ifontdirection) ) then kfontdirection = ifontdirection else kfontdirection = 0 end if call EnterCriticalSection( loc(lpCriticalSection) ) iretb = SetBkMode(wnd%hDC, TRANSPARENT) hFont = CreateFont( kfontsize , 10 , kfontdirection , 0 ,FW_DONTCARE , FALSE , FALSE , FALSE , & ANSI_CHARSET , OUT_DEFAULT_PRECIS , & CLIP_DEFAULT_PRECIS , PROOF_QUALITY , & ior(FIXED_PITCH, FF_ROMAN) , NULL ) iretb = SelectObject(wnd%hdc , hFont) iretb = TextOut(wnd%hDC, ix, iy, txt, len_trim(txt)) iretb = SelectObject(wnd%hdc , GetStockObject(SYSTEM_FONT)) iretb = DeleteObject(hFont) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_text !------------------------------------------------------------------------------------- integer function irgb(ir, ig, ib) implicit none integer, intent(in) :: ir, ig, ib irgb = ir + (ig + (ib * 256)) * 256 return end function irgb !Relative coordinate !------------------------------------------------------------------------------------- subroutine gr_axis(xmin, xmax, ymin, ymax) real, intent(in) :: xmin, xmax, ymin, ymax wnd%xmin = xmin wnd%xmax = xmax wnd%ymin = ymin wnd%ymax = ymax return end subroutine gr_axis !---------------------------------------------------------------- subroutine to_ixy(x, y, ix, iy) real , intent(in ) :: x, y integer, intent(out) :: ix, iy ix = int( wnd%nsize_x * (x - wnd%xmin) / (wnd%xmax - wnd%xmin) ) iy = -int( wnd%nsize_y * (y - wnd%ymin) / (wnd%ymax - wnd%ymin) ) + wnd%nsize_y return end subroutine to_ixy !---------------------------------------------------------------- subroutine gr_wdot(x, y, icol) real, intent(in) :: x, y integer, intent(in) :: icol integer :: ix, iy integer (BOOL):: iretb call to_ixy(x, y, ix, iy) call EnterCriticalSection( loc(lpCriticalSection) ) iretb = SetPixel(wnd%hDC, ix, iy, icol) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_wdot !---------------------------------------------------------------- subroutine gr_wmove(x, y) real, intent(in) :: x, y integer (BOOL):: iretb integer :: ix, iy call to_ixy(x, y, ix, iy) call EnterCriticalSection( loc(lpCriticalSection) ) iretb = MoveToEx(wnd%hDC, ix, iy, NULL) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_wmove !---------------------------------------------------------------- subroutine gr_wline(x, y) real, intent(in) :: x, y integer (BOOL):: iretb integer :: ix, iy call to_ixy(x, y, ix, iy) call EnterCriticalSection( loc(lpCriticalSection) ) iretb = LineTo(wnd%hDC, ix, iy) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_wline !---------------------------------------------------------------- subroutine gr_wtext(x, y, text, icol, ifontsize, ifontdirection) real, intent(in) :: x, y character, intent(in) :: text*(*) integer, optional, intent(in) :: icol integer :: ix, iy integer (BOOL):: iretb integer, optional, intent(in) :: ifontsize, ifontdirection integer :: kfontsize, kfontdirection if ( present(icol) ) iretb = SetTextColor(wnd%hDC, icol) if ( present(ifontsize) ) then kfontsize = ifontsize else kfontsize = 10 end if if ( present(ifontdirection) ) then kfontdirection = ifontdirection else kfontdirection = 0 end if call to_ixy(x, y, ix, iy) call EnterCriticalSection( loc(lpCriticalSection) ) call gr_text(ix, iy, text, icol, kfontsize, kfontdirection) call LeaveCriticalSection( loc(lpCriticalSection) ) return end subroutine gr_wtext !---------------------------------------------------------------- end module uhoplot !===================================================================================== module plotter use uhoplot real :: xpen = 0.0, ypen = 0.0 contains !---------------------------------------------------------------- subroutine move(x, y) real, intent(in) :: x, y xpen = x ypen = y return end subroutine move !---------------------------------------------------------------- subroutine move_rel(dx, dy) real, intent(in) :: dx, dy xpen = xpen + dx ypen = ypen + dy return end subroutine move_rel !---------------------------------------------------------------- subroutine draw(x, y) real, intent(in) :: x, y call gr_wmove(xpen, ypen) call gr_wline(x , y ) xpen = x ypen = y call gr_show() return end subroutine draw !---------------------------------------------------------------- subroutine draw_rel(dx, dy) real, intent(in) :: dx, dy call gr_wmove(xpen, ypen) call gr_wline(xpen + dx, ypen + dy) xpen = xpen + dx ypen = ypen + dy call gr_show() return end subroutine draw_rel !---------------------------------------------------------------- end MODULE plotter