Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/fruit.f90.svn-base

+!------------------------
+! FORTRAN unit test utility
+!
+! Author: Andrew H. Chen meihome @at@ gmail.com
+!------------------------
+!
+! Unit test framework for FORTRAN.  (FoRtran UnIT)
+!
+! This package is to perform unit test for FORTRAN subroutines
+!
+! The method used most are: assert_true, assert_equals
+! 
+! Coding convention:  
+!   1) All methods must be exposed by interface.  i.e. interface init_fruit
+!   2) Variable and methods are lower case connected with underscores.  i.e. init_fruit, and
+!      failed_assert_count
+!
+module fruit
+  use fruit_util
+  implicit none
+  private
+
+  integer, parameter :: MSG_LENGTH = 256
+  integer, parameter :: MAX_MSG_STACK_SIZE = 2000
+  integer, parameter :: MSG_ARRAY_INCREMENT = 50
+  
+  character(*), parameter :: DEFAULT_UNIT_NAME = '_not_set_'
+
+  integer, private, save :: current_max = 50
+  integer, private, save :: successful_assert_count = 0
+  integer, private, save :: failed_assert_count = 0
+  character (len = MSG_LENGTH), private, allocatable :: message_array(:)
+  character (len = MSG_LENGTH), private, save :: msg = '[unit name not set from set_name]: '
+  character (len = MSG_LENGTH), private, save :: unit_name  = DEFAULT_UNIT_NAME
+  integer, private, save :: messageIndex = 1
+
+  integer, private, save :: successful_case_count = 0
+  integer, private, save :: failed_case_count = 0
+  integer, private, save :: testCaseIndex = 1
+  logical, private, save :: last_passed = .false.
+
+  public :: &
+    init_fruit, initializeFruit, fruit_summary, getTestSummary, get_last_message, &
+    is_last_passed, assert_true, assertTrue, assert_equals, assertEquals, &
+    assert_not_equals, assertNotEquals, add_success, addSuccess, &
+    addFail, add_fail, set_unit_name, get_unit_name, &
+    failed_assert_action, get_total_count, getTotalCount, &
+    get_failed_count, getFailedCount, is_all_successful, isAllSuccessful, &
+    run_test_case, runTestCase
+
+  interface initializeFruit
+     module procedure obsolete_initializeFruit_
+  end interface
+
+  interface getTestSummary
+     module procedure obsolete_getTestSummary_  
+  end interface
+
+  interface assertTrue
+     module procedure obsolete_assert_true_logical_
+  end interface
+
+  interface assert_equals
+     module procedure assert_eq_int_
+     module procedure assert_eq_double_
+     module procedure assert_eq_real_
+     module procedure assert_eq_logical_
+     module procedure assert_eq_string_
+     module procedure assert_eq_complex_
+     module procedure assert_eq_real_in_range_
+     module procedure assert_eq_double_in_range_
+
+     module procedure assert_eq_1d_int_
+     module procedure assert_eq_1d_double_
+     module procedure assert_eq_1d_real_
+     module procedure assert_eq_1d_string_
+     module procedure assert_eq_1d_complex_
+     module procedure assert_eq_1d_real_in_range_
+     module procedure assert_eq_1d_double_in_range_
+
+     module procedure assert_eq_2d_int_
+     module procedure assert_eq_2d_double_
+     module procedure assert_eq_2d_real_
+     module procedure assert_eq_2d_complex_
+     module procedure assert_eq_2d_real_in_range_
+     module procedure assert_eq_2d_double_in_range_
+  end interface
+
+  interface assertEquals
+     module procedure assert_eq_int_
+     module procedure assert_eq_double_
+     module procedure assert_eq_real_
+     module procedure assert_eq_logical_
+     module procedure assert_eq_string_
+     module procedure assert_eq_complex_
+     module procedure assert_eq_real_in_range_
+     module procedure assert_eq_double_in_range_
+
+     module procedure assert_eq_1d_int_
+     module procedure assert_eq_1d_double_
+     module procedure assert_eq_1d_real_
+     module procedure assert_eq_1d_string_
+     module procedure assert_eq_1d_complex_
+     module procedure assert_eq_1d_real_in_range_
+     module procedure assert_eq_1d_double_in_range_
+
+     module procedure assert_eq_2d_int_
+     module procedure assert_eq_2d_double_
+     module procedure assert_eq_2d_real_
+     module procedure assert_eq_2d_complex_
+     module procedure assert_eq_2d_real_in_range_
+     module procedure assert_eq_2d_double_in_range_
+  end interface
+
+  interface assert_not_equals
+     module procedure assert_not_equals_real_
+     module procedure assert_not_equals_1d_real_
+     module procedure assert_not_equals_double_
+  end interface
+
+  interface assertNotEquals
+     module procedure assert_not_equals_real_
+     module procedure assert_not_equals_1d_real_
+     module procedure assert_not_equals_double_
+  end interface
+
+  interface addSuccess
+     module procedure obsolete_addSuccess_
+  end interface
+
+  interface add_fail
+     module procedure add_fail_
+     module procedure add_fail_unit_
+  end interface
+
+  interface addFail
+     module procedure add_fail_
+     module procedure add_fail_unit_
+  end interface
+
+  interface getTotalCount
+     module procedure obsolete_getTotalCount_
+  end interface
+
+  interface getFailedCount
+     module procedure obsolete_getFailedCount_
+  end interface
+
+  interface isAllSuccessful
+     module procedure obsolete_isAllSuccessful_
+  end interface
+
+  interface run_test_case
+     module procedure run_test_case_
+     module procedure run_test_case_named_
+  end interface
+
+  interface runTestCase
+     module procedure run_test_case_
+     module procedure run_test_case_named_
+  end interface
+
+contains
+
+  subroutine init_fruit
+    successful_assert_count = 0
+    failed_assert_count = 0
+    messageIndex = 1
+    write (*,*)
+    write (*,*) "Test module initialized"
+    write (*,*)
+    write (*,*) "   . : successful assert,   F : failed assert "
+    write (*,*)
+    if ( .not. allocated(message_array) ) then
+      allocate(message_array(MSG_ARRAY_INCREMENT)) 
+    end if
+  end subroutine init_fruit
+
+  subroutine obsolete_initializeFruit_
+    call obsolete_ ("initializeFruit is OBSOLETE.  replaced by init_fruit")
+    call init_fruit
+  end subroutine obsolete_initializeFruit_
+
+  subroutine obsolete_getTestSummary_
+    call obsolete_ ( "getTestSummary is OBSOLETE.  replaced by fruit_summary")
+    call fruit_summary
+  end subroutine obsolete_getTestSummary_
+
+  ! Run a named test case
+  subroutine run_test_case_named_( tc, tc_name )
+    interface
+       subroutine tc()
+       end subroutine
+    end interface
+    character(*), intent(in) :: tc_name
+
+    integer :: initial_failed_assert_count
+
+    initial_failed_assert_count = failed_assert_count
+
+    ! Set the name of the unit test
+    call set_unit_name( tc_name )
+
+    last_passed = .true.
+
+    call tc()
+
+    if ( initial_failed_assert_count .eq. failed_assert_count ) then
+       ! If no additional assertions failed during the run of this test case
+       ! then the test case was successful
+       successful_case_count = successful_case_count+1
+    else
+       failed_case_count = failed_case_count+1
+    end if
+
+    testCaseIndex = testCaseIndex+1
+    
+    ! Reset the name of the unit test back to the default
+    call set_unit_name( DEFAULT_UNIT_NAME )
+
+  end subroutine run_test_case_named_
+
+  ! Run an 'unnamed' test case
+  subroutine run_test_case_( tc )
+    interface
+       subroutine tc()
+       end subroutine
+    end interface
+
+    call run_test_case_named_( tc, '_unnamed_' )
+
+  end subroutine run_test_case_
+
+  subroutine fruit_summary
+    integer :: i
+
+    write (*,*)
+    write (*,*)
+    write (*,*) '    Start of FRUIT summary: '
+    write (*,*)
+
+    if (failed_assert_count > 0) then
+       write (*,*) 'Some tests failed!'
+    else
+       write (*,*) 'SUCCESSFUL!'
+    end if
+
+    write (*,*)
+    if ( messageIndex > 1) then
+       write (*,*) '  -- Failed assertion messages:'
+
+       do i = 1, messageIndex - 1
+          write (*,"(A)") '   '//trim(strip(message_array(i)))
+       end do
+
+       write (*,*) '  -- end of failed assertion messages.'
+       write (*,*)
+    else
+       write (*,*) '  No messages '
+    end if
+
+    if (successful_assert_count + failed_assert_count /= 0) then
+
+       write (*,*) 'Total asserts :   ', successful_assert_count + failed_assert_count
+       write (*,*) 'Successful    :   ', successful_assert_count
+       write (*,*) 'Failed        :   ', failed_assert_count
+       write (*,'("Successful rate:   ",f6.2,"%")')  real(successful_assert_count) * 100.0 / &
+            real (successful_assert_count + failed_assert_count)
+       write (*, *)
+       write (*,*) 'Successful asserts / total asserts : [ ',&
+            successful_assert_count, '/', successful_assert_count + failed_assert_count, ' ]'
+       write (*,*) 'Successful cases   / total cases   : [ ', successful_case_count, '/', &
+            successful_case_count + failed_case_count, ' ]'
+       write (*, *) '  -- end of FRUIT summary'
+
+    end if
+  end subroutine fruit_summary
+
+  subroutine obsolete_addSuccess_
+    call obsolete_ ("addSuccess is OBSOLETE.  replaced by add_success")
+    call add_success
+  end subroutine obsolete_addSuccess_
+
+  subroutine add_fail_ (message)
+    character (*), intent (in), optional :: message
+    call failed_assert_action('none', 'none', message)
+  end subroutine add_fail_
+
+  subroutine add_fail_unit_ (unitName, message)
+    character (*), intent (in) :: unitName
+    character (*), intent (in) :: message
+
+    call add_fail_ ("[in " //  unitName // "(fail)]: " //  message)
+  end subroutine add_fail_unit_
+
+  subroutine obsolete_isAllSuccessful_(result)
+    logical, intent(out) :: result
+    call obsolete_ ('subroutine isAllSuccessful is changed to function is_all_successful.')
+    result = (failed_assert_count .eq. 0 )
+  end subroutine obsolete_isAllSuccessful_
+
+  subroutine is_all_successful(result)
+    logical, intent(out) :: result
+    result= (failed_assert_count .eq. 0 )
+  end subroutine is_all_successful
+
+  subroutine success_mark_
+    write(*,"(A1)",ADVANCE='NO') '.'
+  end subroutine success_mark_
+
+  subroutine failed_mark_
+    write(*,"(A1)",ADVANCE='NO') 'F'
+  end subroutine failed_mark_
+
+  subroutine increase_message_stack_
+    character(len=MSG_LENGTH) :: msg_swap_holder(current_max)
+
+    if (messageIndex > MAX_MSG_STACK_SIZE) then
+       write(*,*) "Stop because there are too many error messages to put into stack."
+       write (*,*) "Try to increase MAX_MSG_STACK_SIZE if you really need so."
+       call getTestSummary ()
+       stop 1
+    end if
+
+    if (messageIndex > current_max) then
+      msg_swap_holder(1:current_max) = message_array(1:current_max)
+      deallocate(message_array)
+      current_max = current_max + MSG_ARRAY_INCREMENT
+      allocate(message_array(current_max))
+      message_array(1:current_max - MSG_ARRAY_INCREMENT) &
+                   = msg_swap_holder(1: current_max - MSG_ARRAY_INCREMENT)
+    end if
+    
+    message_array (messageIndex) = msg
+    messageIndex = messageIndex + 1
+  end subroutine increase_message_stack_
+
+  function get_last_message()
+    character(len=MSG_LENGTH) :: get_last_message
+    if (messageIndex > 1) then
+       get_last_message = strip(message_array(messageIndex-1))
+    else
+       get_last_message = ''
+    end if
+  end function get_last_message
+
+  subroutine obsolete_getTotalCount_ (count)
+    integer, intent (out) :: count
+    call obsolete_ (' getTotalCount subroutine is replaced by function get_total_count')
+    call get_total_count(count)
+  end subroutine obsolete_getTotalCount_
+
+  subroutine get_total_count(count) 
+    integer, intent(out) :: count
+
+    count = successful_assert_count + failed_assert_count
+  end subroutine get_total_count
+
+  subroutine obsolete_getFailedCount_ (count)
+    integer, intent (out) :: count
+
+    call obsolete_ (' getFailedCount subroutine is replaced by function get_failed_count')
+    call get_failed_count (count)
+
+  end subroutine obsolete_getFailedCount_
+
+  subroutine get_failed_count (count)
+    integer, intent(out) :: count
+    count = failed_assert_count
+  end subroutine get_failed_count
+
+  subroutine obsolete_ (message)
+    character (*), intent (in), optional :: message
+    write (*,*) 
+    write (*,*) "<<<<<<<<<<<<<<<<<<<<<<<<<< WARNING from FRUIT >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
+    write (*,*) message
+    write (*,*) 
+    write (*,*) " old calls will be replaced in the next release in Jan 2009"
+    write (*,*) " Naming convention for all the method calls are changed to: first_name from"
+    write (*,*) " firstName.  Subroutines that will be deleted: assertEquals, assertNotEquals,"
+    write (*,*) " assertTrue, addSuccessful, addFail, etc."
+    write (*,*) "<<<<<<<<<<<<<<<<<<<<<<<<<< WARNING from FRUIT >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
+    write (*,*) 
+  end subroutine obsolete_
+
+  subroutine add_success
+    successful_assert_count = successful_assert_count + 1
+    last_passed = .true.
+    call success_mark_  
+  end subroutine add_success
+
+  subroutine failed_assert_action (expected, got, message)
+    character(*), intent(in) :: expected, got
+    character(*), intent(in), optional :: message
+
+    call make_error_msg_ (expected, got, message)
+    call increase_message_stack_
+    failed_assert_count = failed_assert_count + 1
+    last_passed = .false.
+    call failed_mark_
+  end subroutine failed_assert_action
+
+  subroutine set_unit_name(value)
+    character(*), intent(in) :: value
+    unit_name = strip(value)
+  end subroutine set_unit_name
+
+  subroutine get_unit_name(value)
+    character(*), intent(out) :: value
+    value = strip(unit_name)
+  end subroutine get_unit_name
+
+  subroutine make_error_msg_ (var1, var2, message)
+    character(*), intent(in) :: var1, var2
+    character(*), intent(in), optional :: message
+    msg = '[' // trim(strip(unit_name)) // ']: Expected [' // trim(strip(var1)) &
+          // '], Got [' // trim(strip(var2)) // ']'
+    if (present(message)) then
+       msg = msg // '; User message: [' // message // ']'
+    endif
+  end subroutine make_error_msg_
+
+  function is_last_passed()
+    logical:: is_last_passed
+    is_last_passed = last_passed 
+  end function is_last_passed
+
+  !--------------------------------------------------------------------------------
+  ! all assertions
+  !--------------------------------------------------------------------------------
+  subroutine obsolete_assert_true_logical_(var1, message)
+    logical, intent (in) :: var1
+    character (*), intent (in), optional :: message
+
+    call obsolete_ ('assertTrue subroutine is replaced by function assert_true')
+    call assert_true(var1, message)
+  end subroutine obsolete_assert_true_logical_
+
+  subroutine assert_true (var1, message)
+    logical, intent (in) :: var1
+    character (*), intent (in), optional :: message
+
+    if ( var1 .eqv. .true.) then
+       call add_success
+    else
+       call failed_assert_action(to_s(.true.), to_s(var1), message)
+    end if
+  end subroutine assert_true
+
+  subroutine assert_eq_int_ (var1, var2, message)
+    integer, intent(in) :: var1, var2
+    character (*), intent(in), optional :: message
+
+    if ( var1 .eq. var2) then
+       call add_success
+    else
+       call failed_assert_action (to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_eq_int_
+
+  subroutine assert_eq_logical_ (var1, var2, message)
+    logical, intent (in)  :: var1, var2
+    character (*), intent (in), optional :: message
+
+    if ( var1 .eqv. var2 ) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_eq_logical_
+
+  subroutine assert_eq_string_ (var1, var2, message)
+    character(*), intent (in)  :: var1, var2
+    character (*), intent (in), optional :: message
+
+    if ( trim(strip(var1)) == trim(strip(var2))) then
+       call add_success
+    else
+       call failed_assert_action(var1, var2, message)
+    end if
+  end subroutine assert_eq_string_
+
+  subroutine assert_eq_real_ (var1, var2, message)
+    real, intent (in) :: var1, var2
+    character (*), intent (in), optional :: message
+
+    if ( var1 .eq. var2) then
+       call add_success
+    else
+7      call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_eq_real_
+
+  subroutine assert_eq_double_ (var1, var2, message)
+    double precision, intent (in) :: var1, var2
+    character(*), intent(in), optional :: message
+
+    if ( var1 .eq. var2) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_eq_double_
+
+  subroutine assert_eq_complex_ (var1, var2, message)
+    complex(kind=kind(1.0D0)), intent(IN) :: var1, var2
+    character (*),             intent(IN), optional :: message
+    integer count
+
+    if ( var1 .ne. var2) then
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    else
+       call add_success
+    end if
+
+  end subroutine assert_eq_complex_
+
+  subroutine assert_eq_real_in_range_(var1, var2, var3, message)
+    real, intent (in) :: var1, var2, var3
+    character(*), intent(in), optional :: message
+
+    if ( abs( var1 - var2) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+
+  end subroutine assert_eq_real_in_range_
+
+  subroutine assert_eq_double_in_range_(var1, var2, var3, message)
+    double precision, intent (in) :: var1, var2, var3
+    character(*), intent(in), optional :: message
+
+    if ( abs( var1 - var2) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_eq_double_in_range_
+
+  subroutine assert_eq_1d_int_ (var1, var2, n, message)
+    integer, intent (in) :: n
+    integer, intent (in) :: var1(n), var2(n)
+    character (*), intent (in), optional :: message
+
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)), message)
+          return
+       end if
+    end do loop_dim1
+
+    call add_success
+  end subroutine assert_eq_1d_int_
+
+  subroutine assert_eq_1d_string_ (var1, var2, n, message)
+    integer, intent (in) :: n
+    character(*), intent (in) :: var1(n), var2(n)
+    character (*), intent (in), optional :: message
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( strip(var1(count)) .ne. strip(var2(count))) then
+          call failed_assert_action(var1(count), var2(count), message)
+          return
+       end if
+    end do loop_dim1
+
+    call add_success
+  end subroutine assert_eq_1d_string_
+
+  subroutine assert_eq_1d_real_in_range_(var1, var2, n, var3, message)
+    integer, intent(in) :: n
+    real, intent (in) :: var1(n), var2(n), var3
+    character(*), intent(in), optional :: message
+
+    if ( maxval( abs( var1 - var2)) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1(1)), to_s(var2(1)), &
+                                 '1D array real has difference' // ' ' // message)
+    end if
+  end subroutine assert_eq_1d_real_in_range_
+
+  subroutine assert_eq_1d_double_in_range_(var1, var2, n, var3, message)
+    integer, intent(in) :: n
+    double precision, intent (in) :: var1(n), var2(n), var3
+    character(*), intent(in), optional :: message
+
+    if ( maxval( abs( var1 - var2)) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1(1)), to_s(var2(1)), message)
+    end if
+  end subroutine assert_eq_1d_double_in_range_
+
+  subroutine assert_eq_1d_double (var1, var2, n, message)
+    integer, intent (in) :: n
+    double precision, intent (in) :: var1(n), var2(n)
+    character(*), intent(in), optional :: message
+
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)), &
+               'Array different at count: ' // to_s(count) // ' ' // message)
+          return
+       end if
+    end do loop_dim1
+
+    call add_success
+  end subroutine assert_eq_1d_double
+
+  subroutine assert_eq_2d_real (var1, var2, n, m)
+    integer, intent (in) :: n, m
+    real, intent (in) :: var1(n,m), var2(n,m)
+
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), to_s(var2(count1, count2)),&
+                  'Array (' // to_s(count1) // ',' // to_s( count2) //')')
+             return
+          end if
+       end do loop_dim1
+    end do loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_real
+  
+  subroutine assert_eq_2d_real_in_range_ (var1, var2, n, m, var3, message)
+    integer, intent (in) :: n,m
+    real, intent (in) :: var1(n,m), var2(n,m), var3
+    character(*), intent(in), optional :: message
+
+    if ( maxval( abs( var1 - var2)) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1(1,1)), to_s(var2(1,1)), &
+                                 '2D array real has difference' // ' ' // message)
+    end if
+  
+  end subroutine assert_eq_2d_real_in_range_
+
+  subroutine assert_eq_2d_double (var1, var2, n, m)
+    integer, intent (in) :: n, m
+    double precision, intent (in) :: var1(n,m), var2(n,m)
+
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), to_s(var2(count1, count2)), &
+                  'Array difference at (' // to_s(count1) // ',' // to_s(count2) // ')')
+             return
+          end if
+       end do loop_dim1
+    end do loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_double
+  
+  subroutine assert_eq_2d_double_in_range_ (var1, var2, n, m, var3, message)
+    integer, intent (in) :: n,m
+    double precision, intent (in) :: var1(n,m), var2(n,m), var3
+    character(*), intent(in), optional :: message
+
+    if ( maxval( abs( var1 - var2)) .le. var3) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1(1,1)), to_s(var2(1,1)), &
+                                 '2D array real has difference' // ' ' // message)
+    end if
+  
+  end subroutine assert_eq_2d_double_in_range_
+
+  subroutine assert_eq_2d_int_ (var1, var2, n, m, message)
+    integer, intent (in) :: n, m
+    integer, intent (in) :: var1(n,m), var2(n,m)
+    character (*), intent (in), optional :: message
+
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), &
+                                       to_s(var2(count1, count2)), message)
+             return
+          end if
+       end do loop_dim1
+    end do loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_int_
+
+  subroutine assert_eq_1d_real_ (var1, var2, n, message)
+    integer, intent (in) :: n
+    real, intent (in) :: var1(n), var2(n)
+    character (*), intent (in), optional :: message
+
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)), message)
+          return
+       end if
+    end do loop_dim1
+    call add_success
+  end subroutine assert_eq_1d_real_
+
+  subroutine assert_eq_2d_real_ (var1, var2, n, m, message)
+    integer, intent (in) :: n, m
+    real, intent (in) :: var1(n,m), var2(n,m)
+    character (*), intent(in), optional :: message
+
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), &
+                                       to_s(var2(count1, count2)), message)
+             return
+          end if
+       end do loop_dim1
+    end do loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_real_
+
+  subroutine assert_eq_1d_double_ (var1, var2, n, message)
+    integer, intent (in) :: n
+    double precision, intent (in) :: var1(n), var2(n)
+    character (*), intent (in), optional :: message
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)), message)
+          return
+       end if
+    end do loop_dim1
+
+    call add_success
+  end subroutine assert_eq_1d_double_
+
+  subroutine assert_eq_2d_double_ (var1, var2, n, m, message)
+    integer, intent (in) :: n, m
+    double precision, intent (in) :: var1(n,m), var2(n,m)
+    character (*), intent (in), optional :: message
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), &
+                                       to_s(var2(count1, count2)), message)
+             return
+          end if
+       end do loop_dim1
+    end do loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_double_
+
+  subroutine assert_eq_1d_complex_ (var1, var2, n, message)
+    integer,                   intent(IN) :: n
+    complex(kind=kind(1.0D0)), intent(IN) :: var1(n), var2(n)
+    character (*),             intent(IN), optional :: message
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)), message)
+          return
+       end if
+    enddo loop_dim1
+
+    call add_success
+  end subroutine assert_eq_1d_complex_
+
+  subroutine assert_eq_2d_complex_ (var1, var2, n, m, message)
+    integer,                   intent(IN) :: n, m
+    complex(kind=kind(1.0D0)), intent(IN) :: var1(n,m), var2(n,m)
+    character (*),             intent(IN), optional :: message
+    integer count1, count2
+
+    loop_dim2: do count2 = 1, m
+       loop_dim1: do count1 = 1, n
+          if ( var1(count1,count2) .ne. var2(count1,count2)) then
+             call failed_assert_action(to_s(var1(count1, count2)), &
+                                       to_s(var2(count1, count2)), message)
+             return
+          endif
+       enddo loop_dim1
+    enddo loop_dim2
+
+    call add_success
+  end subroutine assert_eq_2d_complex_
+
+  subroutine assert_not_equals_real_ (var1, var2, message)
+    real, intent (in) :: var1, var2
+    character (*), intent (in), optional :: message
+
+    if ( var1 .ne. var2) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_not_equals_real_
+
+  subroutine assert_not_equals_double_ (var1, var2, message)
+    double precision, intent (in) :: var1, var2
+    character(*), intent(in), optional :: message
+
+    if ( var1 .ne. var2) then
+       call add_success
+    else
+       call failed_assert_action(to_s(var1), to_s(var2), message)
+    end if
+  end subroutine assert_not_equals_double_
+
+  subroutine assert_not_equals_1d_real_ (var1, var2, n)
+    integer, intent (in) :: n
+    real, intent (in) :: var1(n), var2(n)
+
+    integer count
+
+    loop_dim1: do count = 1, n
+       if ( var1(count) .ne. var2(count)) then
+          call failed_assert_action(to_s(var1(count)), to_s(var2(count)),&
+               'Array (' // to_s(count)//')')
+          return
+       end if
+    end do loop_dim1
+
+    call add_success
+
+  end subroutine assert_not_equals_1d_real_
+
+end module fruit

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/fruit_util.f90.svn-base

+module fruit_util
+  private
+  
+  public :: equals, to_s, strip
+  
+  interface equals
+     module procedure equalEpsilon
+     module procedure floatEqual
+     module procedure integerEqual
+     module procedure doublePrecisionEqual
+     module procedure stringEqual
+     module procedure logicalEqual
+  end interface
+
+  interface to_s
+     module procedure to_s_int_
+     module procedure to_s_real_
+     module procedure to_s_logical_
+     module procedure to_s_double_
+     module procedure to_s_complex_
+     module procedure to_s_double_complex_
+     module procedure to_s_string_
+  end interface
+
+contains
+
+  function to_s_int_ (value)
+    implicit none
+    character(len=500):: to_s_int_
+    integer, intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_int_ = adjustl(trim(result))
+  end function to_s_int_
+
+  function to_s_real_ (value)
+    implicit none
+    character(len=500):: to_s_real_
+    real, intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_real_ = adjustl(trim(result))
+  end function to_s_real_
+
+  function to_s_double_ (value)
+    implicit none
+    character(len=500):: to_s_double_
+    double precision, intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_double_ = adjustl(trim(result))
+  end function to_s_double_
+
+  function to_s_complex_ (value)
+    implicit none
+    character(len=500):: to_s_complex_
+    complex, intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_complex_ = adjustl(trim(result))
+  end function to_s_complex_
+
+  function to_s_double_complex_ (value)
+    implicit none
+    character(len=500):: to_s_double_complex_
+    complex(kind=kind(1.0D0)), intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_double_complex_ = adjustl(trim(result))
+  end function to_s_double_complex_
+
+  function to_s_logical_ (value)
+    implicit none
+    character(len=500):: to_s_logical_
+    logical, intent(in) :: value
+    character(len=500) :: result
+    write (result, *) value
+    to_s_logical_ = adjustl(trim(result))
+  end function to_s_logical_
+
+  function to_s_string_ (value)
+    implicit none
+    character(len=500):: to_s_string_
+    character(len=*), intent(in) :: value
+    to_s_string_ = value
+  end function to_s_string_
+
+  function strip(value)
+    implicit none
+    character(len=500):: strip
+    character(len=*), intent(in) :: value
+    strip = trim(adjustl(value))
+  end function strip
+
+  !------------------------
+  ! test if 2 values are close
+  !------------------------
+  !logical function equals (number1, number2) 
+  !  real,  intent (in) :: number1, number2
+  !  
+  !  return equalEpsilon (number1, number2, epsilon(number1))
+  !
+  !end function equals
+
+
+  function equalEpsilon (number1, number2, epsilon ) result (resultValue)
+    real , intent (in) :: number1, number2, epsilon 
+    logical :: resultValue 
+
+    resultValue = .false.
+
+    ! test very small number1
+    if ( abs(number1) < epsilon .and.  abs(number1 - number2) < epsilon ) then
+       resultValue = .true.
+    else 
+       if ((abs(( number1 - number2)) / number1) < epsilon ) then
+          resultValue = .true.
+       else
+          resultValue = .false.
+       end if
+    end if
+
+  end function equalEpsilon
+
+  function floatEqual (number1, number2 ) result (resultValue)
+    real , intent (in) :: number1, number2
+    real :: epsilon 
+    logical :: resultValue 
+
+    resultValue = .false.
+    epsilon = 1E-6
+
+    ! test very small number1
+    if ( abs(number1) < epsilon .and.  abs(number1 - number2) < epsilon ) then
+       resultValue = .true.
+    else 
+       if ((abs(( number1 - number2)) / number1) < epsilon ) then
+          resultValue = .true.
+       else
+          resultValue = .false.
+       end if
+    end if
+  end function floatEqual
+
+  function doublePrecisionEqual (number1, number2 ) result (resultValue)
+    double precision , intent (in) :: number1, number2
+    real :: epsilon 
+    logical :: resultValue 
+
+    resultValue = .false.
+    epsilon = 1E-6
+    !epsilon = epsilon (number1)
+
+    ! test very small number1
+    if ( abs(number1) < epsilon .and.  abs(number1 - number2) < epsilon ) then
+       resultValue = .true.
+    else 
+       if ((abs(( number1 - number2)) / number1) < epsilon ) then
+          resultValue = .true.
+       else
+          resultValue = .false.
+       end if
+    end if
+  end function doublePrecisionEqual
+
+  function integerEqual (number1, number2 ) result (resultValue)
+    integer , intent (in) :: number1, number2
+    logical :: resultValue 
+
+    resultValue = .false.
+
+    if ( number1 .eq. number2 ) then
+       resultValue = .true.
+    else 
+       resultValue = .false.
+    end if
+  end function integerEqual
+
+  function stringEqual (str1, str2 ) result (resultValue)
+    character(*) , intent (in) :: str1, str2
+    logical :: resultValue 
+
+    resultValue = .false.
+
+    if ( str1 .eq. str2 ) then
+       resultValue = .true.
+    end if
+  end function stringEqual
+
+  function logicalEqual (l1, l2 ) result (resultValue)
+    logical, intent (in) :: l1, l2
+    logical              :: resultValue 
+
+    resultValue = .false.
+
+    if ( l1 .eqv. l2 ) then
+       resultValue = .true.
+    end if
+  end function logicalEqual
+
+end module fruit_util

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/integration.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for Body Element Integration
+! ------------------------------------------------------------------------
+MODULE body_element_integration
+USE pde_solver_constants          !> Module for PDE Solver constants (for Gaussian quadrature)
+USE constitutive                  !> Constitutive Matrix module
+USE body_element_interpolation    !> Body Element Interpolation module (shape functions)
+USE dense_matrix_def              !> Dense Matrix data type
+USE field_data                    !> Field Data module
+USE kinematic                     !> Kinematic Matrix module
+USE material_data                 !> Material Property Data module
+USE vector_def                    !> Vector data type
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: bint_emass, bint_estiff, bint_eacc, bint_estress, bint_estrain
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to element mass matrix
+INTERFACE bint_emass
+  MODULE PROCEDURE bint_emass_3_noded_triangular_
+END INTERFACE bint_emass
+
+!> \brief Interface to element stiffness matrix
+INTERFACE bint_estiff
+  MODULE PROCEDURE bint_estiff_3_noded_triangular_
+END INTERFACE bint_estiff
+
+!> \brief Interface to element force vector due to body acceleration
+INTERFACE bint_eacc
+  MODULE PROCEDURE bint_eacc_3_noded_triangular_
+END INTERFACE bint_eacc
+
+!> \brief Interface to element force vector due to initial stress
+INTERFACE bint_estress
+  MODULE PROCEDURE bint_estress_3_noded_triangular_
+END INTERFACE bint_estress
+
+!> \brief Interface to element force vector due to initial strain
+INTERFACE bint_estrain
+  MODULE PROCEDURE bint_estrain_3_noded_triangular_
+END INTERFACE bint_estrain
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element mass matrix (3-noded triangular)
+!!
+!!  \param  iel   Element number
+!!  \param  emass Element mass matrix
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_A N^T * rho * N dA
+!!
+!!  where N is the shape function matrix, rho is the material density,
+!!  and A is the area of the element. The routine uses the technique of
+!!  Gaussian quadrature to numerically evaluate the integral. However,
+!!  note that the integration points and weights are selected such that
+!!  the exact integral is obtained. Consequently, only the usual floating
+!!  point error is present in the result of the integration.
+! ------------------------------------------------------------------------
+SUBROUTINE bint_emass_3_noded_triangular_ (iel, emass)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(matrixT), INTENT(INOUT) :: emass
+  TYPE(matrixT) :: N        !> shape function matrix
+  DOUBLE PRECISION :: rho   !> material density
+  DOUBLE PRECISION :: A     !> element area
+  INTEGER :: m    !> material type
+  INTEGER :: i    !> loop variable
+
+  m = fld_getMaterial(iel)    !> get material type
+  rho = mtl_getDens(m)        !> get density
+  A = fld_volElem(iel)        !> get element area
+
+  !> initialize element mass matrix
+  CALL dm_init(emass, NNODEL*NDIM,NNODEL*NDIM)
+
+  !> perform Gaussian quadrature
+  DO i = 1,NGAUSS_ELEM
+
+    !> get shape functions at integration point
+    CALL bshp_shape(GAUSS_PT_ELEM(i,1),GAUSS_PT_ELEM(i,2), N)
+
+    !> increment element mass matrix
+    emass%dat = emass%dat &
+                + GAUSS_WT_ELEM(i)*rho*A*MATMUL(TRANSPOSE(N%dat),N%dat)
+
+  END DO
+
+  !> deallocate shape function matrix
+  CALL dm_clean(N)
+
+END SUBROUTINE bint_emass_3_noded_triangular_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element stiffness matrix (3-noded triangular)
+!!
+!!  \param  iel     Element number
+!!  \param  estiff  Element stiffness matrix
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_A B^T * D * B dA
+!!
+!!  where B is the kinematic matrix, D is the constitutive matrix,
+!!  and A is the area of the element. Since B and D are constants for this
+!!  type of element, the integration reduces to (B^T * D * B * A).
+! ------------------------------------------------------------------------
+SUBROUTINE bint_estiff_3_noded_triangular_ (iel, estiff)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(matrixT), INTENT(INOUT) :: estiff
+  TYPE(matrixT) :: B      !> kinematic matrix
+  TYPE(matrixT) :: D      !> constitutive matrix
+  DOUBLE PRECISION :: emod, nu    !> material properties
+  DOUBLE PRECISION :: A           !> element area
+  INTEGER :: m    !> material type
+
+  m = fld_getMaterial(iel)    !> get material type
+  emod = mtl_getEmod(m)       !> get elastic modulus
+  nu = mtl_getPois(m)         !> get Poisson's ratio
+  A = fld_volElem(iel)        !> get element area
+
+  !> initialize element stiffness matrix
+  CALL dm_init(estiff, NNODEL*NDIM,NNODEL*NDIM)
+
+  !> get kinematic and constitutive matrices
+  CALL bmatrix(iel, B)
+  CALL dmatrix(emod,nu, D)
+
+  !> compute element stiffness matrix
+  estiff%dat = A * MATMUL(TRANSPOSE(B%dat),MATMUL(D%dat,B%dat))
+
+  !> deallocate kinematic and constitutive matrices
+  CALL dm_clean(B)
+  CALL dm_clean(D)
+
+END SUBROUTINE bint_estiff_3_noded_triangular_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element force vector due to body acceleration (3-noded triangular)
+!!
+!!  \param  iel     Element number
+!!  \param  eload   Element load vector
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_A N^T * rho * f dA
+!!
+!!  where N is the shape function matrix, rho is the material density,
+!!  f is the element body acceleration vector, and A is the area of the
+!!  element. The routine uses the technique of Gaussian quadrature to
+!!  numerically evaluate the integral. However, note that the integration
+!!  points and weights are selected such that the exact integral is
+!!  obtained. Consequently, only the usual floating point error is present
+!!  in the result of the integration.
+! ------------------------------------------------------------------------
+SUBROUTINE bint_eacc_3_noded_triangular_ (iel, eload)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(vectorT), INTENT(INOUT) :: eload
+  TYPE(matrixT) :: N        !> shape function matrix
+  TYPE(vectorT) :: f        !> body acceleration vector
+  DOUBLE PRECISION :: rho   !> density
+  DOUBLE PRECISION :: A     !> element area
+  INTEGER :: m    !> material type
+  INTEGER :: i,j  !> loop variables
+
+  !> initialize body acceleration vector
+  CALL vec_init(f, NNODEL*NDIM)
+
+  !> load body acceleration at nodes
+  DO i = 1,NNODEL
+    DO j = 1,NDIM
+      CALL vec_set( f, (i-1)*NDIM+j, fld_getBodyAcc( fld_getConnect(iel,i), j ) )
+    END DO
+  END DO
+
+  m = fld_getMaterial(iel)    !> get material type
+  rho = mtl_getDens(m)        !> get density
+  A = fld_volElem(iel)        !> get element area
+
+  !> initialize element load vector
+  CALL vec_init(eload, NNODEL*NDIM)
+
+  !> perform Gaussian quadrature
+  DO i = 1,NGAUSS_ELEM
+
+    !> get shape functions at integration point
+    CALL bshp_shape(GAUSS_PT_ELEM(i,1),GAUSS_PT_ELEM(i,2), N)
+
+    !> increment element load vector
+    eload%dat = eload%dat &
+                + GAUSS_WT_ELEM(i)*rho*A*MATMUL(TRANSPOSE(N%dat),f%dat)
+
+  END DO
+
+  !> deallocate objects
+  CALL vec_clean(f)
+  CALL dm_clean(N)
+
+END SUBROUTINE bint_eacc_3_noded_triangular_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element force vector due to initial stress (3-noded triangular)
+!!
+!!  \param  iel     Element number
+!!  \param  eload   Element load vector
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_A B^T * sig0 dA
+!!
+!!  where B is the kinematic matrix, sig0 is the initial stress vector,
+!!  and A is the area of the element. Since B and sig0 are constant for
+!!  the 3-noded triangular element, the integration reduces to
+!!  (B^T * sig0 * A).
+! ------------------------------------------------------------------------
+SUBROUTINE bint_estress_3_noded_triangular_ (iel, eload)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(vectorT), INTENT(INOUT) :: eload
+  TYPE(matrixT) :: B        !> kinematic matrix
+  TYPE(vectorT) :: sig0     !> initial stress vector
+  DOUBLE PRECISION :: A     !> element area
+
+  !> initialize initial stress vector
+  CALL vec_init(sig0, NTNS)
+
+  !> load initial stress vector
+  CALL vec_set( sig0, 1, fld_getStressElem(iel,s11) )
+  CALL vec_set( sig0, 2, fld_getStressElem(iel,s22) )
+  CALL vec_set( sig0, 3, fld_getStressElem(iel,s12) )
+
+  !> get element area
+  A = fld_volElem(iel)
+
+  !> get kinematic matrix
+  CALL bmatrix(iel,B)
+
+  !> initialize element load vector
+  CALL vec_init(eload, NNODEL*NDIM)
+
+  !> compute element load vector
+  eload%dat = A * MATMUL(TRANSPOSE(B%dat),sig0%dat)
+
+  !> deallocate objects
+  CALL vec_clean(sig0)
+  CALL dm_clean(B)
+
+END SUBROUTINE bint_estress_3_noded_triangular_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element force vector due to initial strain (3-noded triangular)
+!!
+!!  \param  iel     Element number
+!!  \param  eload   Element load vector
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_A B^T * D * eps0 dA
+!!
+!!  where B is the kinematic matrix, D is the constitutive matrix,
+!!  eps0 is the initial strain vector, and A is the area of the element.
+!!  Since B, D, and eps0 are constant for the 3-noded triangular element,
+!!  the integration reduces to (B^T * D * eps0 * A).
+! ------------------------------------------------------------------------
+SUBROUTINE bint_estrain_3_noded_triangular_ (iel, eload)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(vectorT), INTENT(INOUT) :: eload
+  TYPE(matrixT) :: B        !> kinematic matrix
+  TYPE(matrixT) :: D        !> constitutive matrix
+  TYPE(vectorT) :: eps0     !> initial strain vector
+  INTEGER :: m                  !> material number
+  DOUBLE PRECISION :: emod, nu  !> material properties
+  DOUBLE PRECISION :: A         !> element area
+
+  !> initialize initial strain vector
+  CALL vec_init(eps0, NTNS)
+
+  !> load initial strain vector
+  CALL vec_set( eps0, 1, fld_getStrainElem(iel,s11) )
+  CALL vec_set( eps0, 2, fld_getStrainElem(iel,s22) )
+  CALL vec_set( eps0, 3, fld_getStrainElem(iel,s12) )
+
+  !> get element area
+  A = fld_volElem(iel)
+
+  !> get kinematic matrix
+  CALL bmatrix(iel,B)
+
+  !> get material properties
+  m = fld_getMaterial(iel)
+  emod = mtl_getEmod(m)
+  nu = mtl_getPois(m)
+
+  !> get constitutive matrix
+  CALL dmatrix(emod,nu, D)
+
+  !> initialize element load vector
+  CALL vec_init(eload, NNODEL*NDIM)
+
+  !> compute element load vector
+  eload%dat = A * MATMUL( TRANSPOSE(B%dat), MATMUL(D%dat, eps0%dat) )
+
+  !> deallocate objects
+  CALL vec_clean(eps0)
+  CALL dm_clean(B)
+  CALL dm_clean(D)
+
+END SUBROUTINE bint_estrain_3_noded_triangular_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE body_element_integration
+
+
+
+
+! ------------------------------------------------------------------------
+!>  \brief Module for Traction Element Integration
+! ------------------------------------------------------------------------
+MODULE traction_element_integration
+USE pde_solver_constants              !>
+USE traction_element_interpolation    !> traction element shape functions
+USE dense_matrix_def                  !> Dense Matrix data type
+USE vector_def                        !> Vector data type
+USE boundary_data                     !> boundary element data
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: tint_etrac
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to element mass matrix
+INTERFACE tint_etrac
+  MODULE PROCEDURE tint_etrac_2_noded_linear_
+END INTERFACE tint_etrac
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Element load vector for traction element (2-noded linear)
+!!
+!!  \param  iel   Element number
+!!  \param  eload Element load vector
+!!
+!!  This routine performs the following integration:
+!!
+!!      \int_0^1 Nt^T * T^T * Nt * tt * lt ds
+!!
+!!  where Nt is the shape function matrix, T is the transformation matrix,
+!!  tt is the local traction vector, and lt is the length of the element.
+!!  The routine uses the technique of Gaussian quadrature to numerically
+!!  evaluate the integral. However, note that the integration points and
+!!  weights are selected such that the exact integral is obtained.
+!!  Consequently, only the usual floating point error is present in the
+!!  result of the integration.
+! ------------------------------------------------------------------------
+SUBROUTINE tint_etrac_2_noded_linear_ (iel, eload)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(vectorT), INTENT(INOUT) :: eload
+  TYPE(matrixT) :: Nt       !> shape function matrix
+  TYPE(matrixT) :: T        !> transformation matrix
+  TYPE(vectorT) :: tt       !> local traction vector
+  TYPE(surfLoadT) :: t1,t2  !> nodal tractions
+  DOUBLE PRECISION :: lt    !> element length
+  INTEGER :: i    !> loop variable
+
+  lt = bnd_lenBoundElem(iel)    !> get element length
+
+  !> get transformation matrix
+  CALL tshp_transform(iel, T)
+
+  !> get nodal tractions
+  t1 = bnd_getTrac(iel,1)
+  t2 = bnd_getTrac(iel,2)
+
+  !> set up traction vector
+  CALL vec_init(tt, NNODELB*NDIM)
+  CALL vec_set(tt, 1, t1%sig_nt)
+  CALL vec_set(tt, 2, t1%sig_nn)
+  CALL vec_set(tt, 3, t2%sig_nt)
+  CALL vec_set(tt, 4, t2%sig_nn)
+
+  !> initialize element load vector
+  CALL vec_init(eload, NNODELB*NDIM)
+
+  !> perform Gaussian quadrature
+  DO i = 1,NGAUSS_BOUND
+
+    !> get shape functions at integration point
+    CALL tshp_shape(GAUSS_PT_BOUND(i), Nt)
+
+    !> increment element load vector
+    eload%dat = eload%dat &
+                + GAUSS_WT_BOUND(i)*lt &
+                  * MATMUL(MATMUL(TRANSPOSE(Nt%dat),MATMUL(TRANSPOSE(T%dat),Nt%dat)),tt%dat)
+
+  END DO
+
+  !> deallocate objects
+  CALL dm_clean(Nt)
+  CALL dm_clean(T)
+  CALL vec_clean(tt)
+
+END SUBROUTINE tint_etrac_2_noded_linear_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE traction_element_integration

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/interpolation.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for Body Element Interpolation
+! ------------------------------------------------------------------------
+MODULE body_element_interpolation
+USE system_constants    !> Global constants (for size of shape function matrix)
+USE dense_matrix_def    !> Dense Matrix data type
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: bshp_shape
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to shape function access program
+INTERFACE bshp_shape
+  MODULE PROCEDURE bshp_shape_3_noded_tri_
+END INTERFACE bshp_shape
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Shape function matrix for body element
+!!
+!!  \param  l1  Area coordinate 1
+!!  \param  l2  Area coordinate 2
+!!  \param  N   Shape function matrix
+!!
+!!  This routine builds the shape function matrix for a 3-noded triangular
+!!  body element. It has the following form:
+!!
+!!    N = [ l1   0  l2   0  l3   0 ]
+!!        [  0  l1   0  l2   0  l3 ]
+!!
+!!  where l3 = 1-l1-l2 and l1+l2 <= 1 is assumed.
+! ------------------------------------------------------------------------
+SUBROUTINE bshp_shape_3_noded_tri_ (l1,l2, N)
+  DOUBLE PRECISION, INTENT(IN) :: l1,l2
+  TYPE(matrixT), INTENT(INOUT) :: N
+  DOUBLE PRECISION :: l3    !> Area coordinate 3
+
+  !> compute third area coordinate
+  l3 = 1.d0-l1-l2
+
+  !> initialize shape function matrix
+  CALL dm_init(N, NDIM, NDIM*NNODEL)
+
+  !> build shape function matrix
+  CALL dm_set(N, 1,1, l1)
+  CALL dm_set(N, 2,2, l1)
+  CALL dm_set(N, 1,3, l2)
+  CALL dm_set(N, 2,4, l2)
+  CALL dm_set(N, 1,5, l3)
+  CALL dm_set(N, 2,6, l3)
+
+END SUBROUTINE bshp_shape_3_noded_tri_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE body_element_interpolation
+
+
+
+
+! ------------------------------------------------------------------------
+!>  \brief Module for Traction Element Interpolation
+! ------------------------------------------------------------------------
+MODULE traction_element_interpolation
+USE system_constants    !> Global constants (for size of shape function matrix)
+USE dense_matrix_def    !> Dense Matrix data type
+USE field_data          !> Field Data module
+USE boundary_data          !> Boundary Data module
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: tshp_shape, tshp_transform
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to shape function access program
+INTERFACE tshp_shape
+  MODULE PROCEDURE tshp_shape_2_noded_linear_
+END INTERFACE tshp_shape
+
+!> \brief Interface to transformation matrix access program
+INTERFACE tshp_transform
+  MODULE PROCEDURE tshp_transform_2_noded_linear_
+END INTERFACE tshp_transform
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Shape function matrix for traction element
+!!
+!!  \param  s   Local coordinate (0 <= s <= 1)
+!!  \param  Nt  Shape function matrix
+!!
+!!  This routine builds the shape function matrix for a 2-noded traction
+!!  element. It has the following form:
+!!
+!!    Nt = [ 1-s   0    s   0 ]
+!!         [  0   1-s   0   s ]
+! ------------------------------------------------------------------------
+SUBROUTINE tshp_shape_2_noded_linear_ (s, Nt)
+  DOUBLE PRECISION, INTENT(IN) :: s
+  TYPE(matrixT), INTENT(INOUT) :: Nt
+  DOUBLE PRECISION :: one_minus_s   !> for efficiency
+
+  !> compute 1-s
+  one_minus_s = 1.d0-s
+
+  !> initialize shape function matrix
+  CALL dm_init(Nt, NDIM, NDIM*NNODELB)
+
+  !> build shape function matrix
+  CALL dm_set(Nt, 1,1, one_minus_s)
+  CALL dm_set(Nt, 2,2, one_minus_s)
+  CALL dm_set(Nt, 1,3, s)
+  CALL dm_set(Nt, 2,4, s)
+
+END SUBROUTINE tshp_shape_2_noded_linear_
+
+
+
+! ------------------------------------------------------------------------
+!>  \brief Transformation matrix for traction element
+!!
+!!  \param  i   Element number
+!!  \param  T   Transformation matrix
+!!
+!!  This routine builds the transformation matrix for a 2-noded traction
+!!  element. It determines the coordinates of the element by accessing the
+!!  field data for the problem. The matrix has the following form:
+!!
+!!    T = [  cos(theta)   sin(theta) ]
+!!        [ -sin(theta)   cos(theta) ]
+!!
+!!  where theta = atan( (y2-y1) / (x2-x1) ). The points (x1,y1) and
+!!  (x2,y2) are the beginning and end coordinates of the element,
+!!  respectively.
+! ------------------------------------------------------------------------
+SUBROUTINE tshp_transform_2_noded_linear_ (i, T)
+  INTEGER, INTENT(IN) :: i
+  TYPE(matrixT), INTENT(INOUT) :: T
+  DOUBLE PRECISION :: x1,y1, x2,y2    !> element coordinates
+  DOUBLE PRECISION :: length          !> length of element
+  DOUBLE PRECISION :: c,s             !> for cos(theta) and sin(theta)
+
+  x1 = fld_getCoord( bnd_getConnect(i,1), 1 )
+  y1 = fld_getCoord( bnd_getConnect(i,1), 2 )
+  x2 = fld_getCoord( bnd_getConnect(i,2), 1 )
+  y2 = fld_getCoord( bnd_getConnect(i,2), 2 )
+
+  length = bnd_lenBoundElem(i)
+
+  c = (x2-x1)/length
+  s = (y2-y1)/length
+
+  CALL dm_init(T, NDIM,NDIM)
+
+  CALL dm_set(T, 1,1, c)
+  CALL dm_set(T, 2,1, -s)
+  CALL dm_set(T, 1,2, s)
+  CALL dm_set(T, 2,2, c)
+
+END SUBROUTINE tshp_transform_2_noded_linear_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE traction_element_interpolation

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/kinematic.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for Kinematic Matrix
+! ------------------------------------------------------------------------
+MODULE kinematic
+USE system_constants    !> Global constants (for size of matrix)
+USE dense_matrix_def    !> Dense Matrix data type
+USE field_data          !> Field Data module
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: bmatrix
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to kinematic matrix access program
+INTERFACE bmatrix
+  MODULE PROCEDURE bmatrix_3_noded_small_strain_
+END INTERFACE bmatrix
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Kinematic matrix (3-noded, small strain)
+!!
+!!  \param  iel   Element number
+!!  \param  Bmat  Kinematic matrix
+!!
+!!  This routine builds the kinematic matrix for a 3-noded triangular
+!!  element considering small strain. This is sometimes referred to in the
+!!  finite element literature as the B matrix, which is defined as:
+!!
+!!        Bmat = L*N
+!!
+!!  where L is the associated linear differential operator and N is the
+!!  matrix of shape functions. Using the mathematical definition of
+!!  strain:
+!!
+!!        L = [ d/dx     0   ]
+!!            [    0   d/dy  ]
+!!            [ d/dy   d/dx  ]
+!!
+!!  where d/dx and d/dy are partial derivatives with respect to the
+!!  coordinate axes. Combining this with the shape functions based on
+!!  area coordinates for a 3-noded triangular element, the kinematic
+!!  matrix is:
+!!
+!!        Bmat = ( 1/(2*A) ) *  [ b1    0   b2    0   b3    0  ]
+!!                              [  0   c1    0   c2    0   c3  ]
+!!                              [ c1   b1   c2   b2   c3   b3  ]
+!!
+!!  where bi and ci are defined as:
+!!
+!!        bi = yj - yk      for {i,j,k} \in {{1,2,3},{2,3,1},{3,1,2}}
+!!        ci = xk - xj
+!!
+!!  and (xi,yi) are the nodal coordinates of the element.
+! ------------------------------------------------------------------------
+SUBROUTINE bmatrix_3_noded_small_strain_ (iel, Bmat)
+  INTEGER, INTENT(IN) :: iel
+  TYPE(matrixT), INTENT(INOUT) :: Bmat
+  DOUBLE PRECISION, DIMENSION(NNODEL) :: x,y  !> for nodal coords
+  DOUBLE PRECISION, DIMENSION(NNODEL) :: b,c  !> for coefficients
+  DOUBLE PRECISION :: A       !> area of element
+  DOUBLE PRECISION :: coef    !> area coefficient
+  INTEGER :: i    !> loop variable
+
+  !> get area of element
+  A = fld_volElem(iel)
+  coef = 1.d0/(2.d0*A)
+
+  !> get nodal coords
+  DO i = 1,NNODEL
+    x(i) = fld_getCoord( fld_getConnect(iel,i), 1 )
+    y(i) = fld_getCoord( fld_getConnect(iel,i), 2 )
+  END DO
+
+  !> compute coefficients
+  b(1) = y(2)-y(3)
+  b(2) = y(3)-y(1)
+  b(3) = y(1)-y(2)
+  b = coef*b
+
+  c(1) = x(3)-x(2)
+  c(2) = x(1)-x(3)
+  c(3) = x(2)-x(1)
+  c = coef*c
+
+  !> initialize matrix
+  CALL dm_init(Bmat, NTNS,NDIM*NNODEL)
+
+  !> set matrix contents
+  Bmat%dat(1,1:NDIM*NNODEL:2) = b(1:NNODEL)
+  Bmat%dat(2,2:NDIM*NNODEL:2) = c(1:NNODEL)
+  Bmat%dat(3,1:NDIM*NNODEL:2) = c(1:NNODEL)
+  Bmat%dat(3,2:NDIM*NNODEL:2) = b(1:NNODEL)
+
+END SUBROUTINE bmatrix_3_noded_small_strain_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE kinematic

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/linear_solver.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module defining solver for linear systems of equations
+! ------------------------------------------------------------------------
+MODULE linear_solver
+USE log_message_control   !> Printing log/error messages
+USE log_messages          !> Log/error codes
+USE band_sym_matrix_def   !> Banded Symmetric Matrix data type
+USE vector_def            !> Vector data type
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+!> Exported interfaces
+PUBLIC :: lin_solve
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+!> Sender code for LINear SoLVer module
+INTEGER, PARAMETER :: sdr = LINSLV
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!>  \brief  Interface to linear solver routine
+INTERFACE lin_solve
+  MODULE PROCEDURE lin_solve_
+  MODULE PROCEDURE lin_solve_exc_
+END INTERFACE lin_solve
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Linear system of equations solver (non-exception checking)
+!!
+!!  \param A    Reference a banded symmetric matrix of coefficients
+!!  \param b    Reference to a vector of right-hand-sides
+!!  \param x    Reference to the vector for the solution (need not be initialized)
+!!  \param hbw  Half bandwidth of A (including the diagonal)
+!!  \param n    Number of equations (number of rows in A / number of entries in b)
+!!  \param info Error code for LAPACK solver (info.EQ.0 -> success | i.LT.0 -> illegal input | i.GT.0 -> A not pos def)
+!!
+!!  This routine solves the system A*x = b. The matrix A should be
+!!  symmetric and positive definite. If the matrix A has not been
+!!  previously used in the solution of a system, it computes the Cholesky
+!!  decomposition using DPBTRF from the LAPACK library. Once the
+!!  decomposition is available it solves the system using DPBTRS, also
+!!  from the LAPACK library.
+! ------------------------------------------------------------------------
+SUBROUTINE lin_solve_ (A,b,x, hbw,n, info)
+  TYPE(bandSymMatrixT), INTENT(INOUT) :: A
+  TYPE(vectorT), INTENT(IN) :: b
+  TYPE(vectorT), INTENT(OUT) :: x
+  INTEGER, INTENT(IN) :: hbw,n
+  INTEGER, INTENT(OUT) :: info
+  CHARACTER, PARAMETER :: uplo = 'u'  !> indicates upper band storage
+  INTEGER, PARAMETER :: nrhs = 1      !> number of right-hand-sides
+  INTEGER :: kd                       !> number of superdiagonals
+  INTEGER :: ldab,ldb                 !> leading dimensions of A and b
+  DOUBLE PRECISION, DIMENSION(hbw,n) :: work_mat  !> workspace for matrix
+  DOUBLE PRECISION, DIMENSION(n) :: work_vec      !> workspace for vector
+
+  kd = hbw-1    !> number of off-diagonal bands
+  ldab = hbw    !> leading dimension of A in banded storage format
+  ldb = n       !> leading dimension of right-hand-sides
+
+  !> if A has already been decomposed, just solve the system
+  IF (bsm_isDecomposed(A)) THEN
+
+    work_mat = A%decomp   !> set matrix workspace to Cholesky decomposition
+    work_vec = b%dat      !> set vector workspace to right-hand-side
+
+  !> if A is not yet decomposed, decompose and then solve
+  ELSE
+
+    work_mat = A%dat    !> set matrix workspace to A in banded storage
+
+    !> perform Cholesky decomposition, work_mat is the decomposition on exit
+    CALL DPBTRF( uplo, n,kd, work_mat,ldab, info )
+
+    !> if decomposition was successful, set the decomposition in A
+    IF (info.EQ.0) THEN
+      CALL bsm_setDecomp(A,work_mat)
+
+    !> otherwise, return
+    ELSE
+      RETURN
+    END IF
+
+    work_vec = b%dat  !> set vector workspace to right-hand-side
+
+  END IF
+
+
+  !> solve the system, work_vec is the solution on exit
+  CALL DPBTRS( uplo, n,kd,nrhs, work_mat,ldab, work_vec,ldb, info )
+
+  !> if solution was successful, initialize solution vector and set its data to the solution
+  IF (info.EQ.0) THEN
+    CALL vec_init(x,n)
+    x%dat = work_vec
+  END IF
+
+END SUBROUTINE lin_solve_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Linear system of equations solver (exception checking)
+!!
+!!  \param A    Reference a banded symmetric matrix of coefficients
+!!  \param b    Reference to a vector of right-hand-sides
+!!  \param x    Reference to the vector for the solution (need not be initialized)
+!!  \param hbw  Half bandwidth of A (including the diagonal)
+!!  \param n    Number of equations (number of rows in A / number of entries in b)
+!!  \param info Error code for LAPACK solver (info.EQ.0 -> success | i.LT.0 -> illegal input | i.GT.0 -> A not pos def)
+!!  \param exc  Error code
+!!
+!!  This routine solves the system A*x = b. The matrix A should be
+!!  symmetric and positive definite. If the matrix A has not been
+!!  previously used in the solution of a system, it computes the Cholesky
+!!  decomposition using DPBTRF from the LAPACK library. Once the
+!!  decomposition is available it solves the system using DPBTRS, also
+!!  from the LAPACK library. If the number of rows in A does not match
+!!  the number of entries in x, it returns a DIMEN exception.
+! ------------------------------------------------------------------------
+SUBROUTINE lin_solve_exc_ (A,b,x, hbw,n, info, exc)
+  TYPE(bandSymMatrixT), INTENT(INOUT) :: A
+  TYPE(vectorT), INTENT(IN) :: b
+  TYPE(vectorT), INTENT(OUT) :: x
+  INTEGER, INTENT(IN) :: hbw,n
+  INTEGER, INTENT(OUT) :: info
+  INTEGER, INTENT(OUT) :: exc
+
+  !> make sure number of equations in A matches number of right-hand-sides in b
+  !> also double-check that input size characteristics match those of A
+  IF (bsm_numRows(A).NE.vec_length(b) &
+        .OR. bsm_halfBW(A).NE.hbw &
+        .OR. bsm_numRows(A).NE.n ) THEN
+    exc=DIMEN
+    CALL log_printLogMsg(exc,sdr)
+    info=-1
+    RETURN
+  ELSE
+    exc=OK
+  END IF
+
+  !> call the non-exception checking solver
+  CALL lin_solve(A,b,x, hbw,n, info)
+
+  !> if solution was unsuccessful, it means A was not positive definite
+  !! (since the selected LAPACK solver only solves systems with banded
+  !! symmetric positive definite coefficient matrices, raise POSDEF
+  !! exception in this case)
+  IF (info.NE.0) THEN
+    CALL vec_clean(x)   !> destroy x since it is not meaningful
+    exc=POSDEF
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+  ELSE
+    exc=OK
+  END IF
+
+END SUBROUTINE lin_solve_exc_
+
+
+END MODULE linear_solver
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/linear_solver_test.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for testing Linear Solver module
+! ------------------------------------------------------------------------
+MODULE linear_solver_test
+USE fruit                   !> Unit testing framework
+USE log_message_control     !> Printing log/error messages
+USE log_messages            !> Log/error codes
+USE band_sym_matrix_def     !> Banded Symmetric Matrix data type
+USE vector_def              !> Vector data type
+USE linear_solver           !> Linear Solver module
+
+IMPLICIT NONE
+
+
+CONTAINS
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for DIMEN exception
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix
+!!  \param test_b   Test right-hand-side vector
+!!  \param test_x   Dummy vector for solution
+!!  \param testName Filename for log file (required for exceptions)
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expMsg   Expected log message
+!!  \param expInfo  Expected LAPACK info code
+!!  \param actMsg   Actual log message
+!!  \param actInfo  Actual LAPACK info code
+!!
+!!  This test tries to solve a system of equations when the dimensions of
+!!  the coefficient matrix and right-hand-side vector are not compatible.
+!!  That is, the number of rows is not the same. It should throw a DIMEN
+!!  exception in this case.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_DIMEN
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_DIMEN'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: hbw=3,n=10
+  INTEGER, PARAMETER :: expMsg=DIMEN, expInfo=-1
+  INTEGER :: actMsg, actInfo
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize matrix and vector (data contents not important)
+  CALL bsm_init(test_A, hbw,n)
+  CALL vec_init(test_b, n+1)    !> note differing dimension
+
+  !> attempt to solve the system
+  CALL lin_solve(test_A,test_b,test_x, hbw,n, actInfo, actMsg)
+  CALL assertEquals(expInfo,actInfo)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate objects
+  CALL log_closeLogFile()
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_DIMEN
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSDEF exception
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix (not positive definite)
+!!  \param test_b   Test right-hand-side vector
+!!  \param test_x   Dummy vector for solution
+!!  \param testName Filename for log file (required for exceptions)
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expMsg   Expected log message
+!!  \param expInfo  Expected LAPACK info code
+!!  \param actMsg   Actual log message
+!!  \param actInfo  Actual LAPACK info code
+!!
+!!  This test tries to solve a system of equations when the coefficient
+!!  matrix, test_A, is not positive definite. It should throw a POSDEF
+!!  exception in this case.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_POSDEF
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_POSDEF'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: hbw=3,n=10
+  INTEGER, PARAMETER :: expMsg=POSDEF, expInfo=1  !> matrix minor 1 is not pos def
+  INTEGER :: actMsg, actInfo
+  INTEGER :: i,j    !> loop variables
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> build a symmetric, but not positive definite, matrix in test_A
+  !! (Matrix is all -2 on the bands, therefore some of the matrix minors
+  !! do not have determinant greater than 0. In particular, matrix minor
+  !! 1 has negative determinant)
+  CALL bsm_init(test_A, hbw,n)
+  DO i = 1,n
+    DO j = i,MIN(i+hbw-1,n)
+      CALL bsm_set(test_A, i,i, -2.d0)
+    END DO
+  END DO
+
+  !> build test_b = [1..n]
+  CALL vec_init(test_b, n)
+  DO i = 1,n
+    CALL vec_set(test_b, i, DBLE(i))
+  END DO
+
+  !> attempt to compute solution
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, actInfo, actMsg)
+  CALL assertEquals(expInfo,actInfo)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate objects
+  CALL log_closeLogFile()
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_POSDEF
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct solution to linear system of equations (general case)
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix
+!!  \param test_b   Test right-hand-side vector
+!!  \param test_x   Vector for solution
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expDat   Expected solution
+!!  \param toler    Tolerance on solution (floating point error)
+!!
+!!  This test solves a system of equations, test_A*test_x = test_b, for
+!!  test_x in the general case for a banded symmetric positive definite
+!!  matrix, test_A. The right-hand-sides, test_b, are computed through
+!!  multiplication of test_A and test_x. That is, the input is
+!!  test_b = test_A*test_x. Therefore, the expected output from the solver
+!!  with inputs test_A and test_b is the original vector test_x.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_VAL
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_VAL'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  INTEGER, PARAMETER :: hbw=3,n=10
+  DOUBLE PRECISION, DIMENSION(n) :: expDat
+  DOUBLE PRECISION, PARAMETER :: toler = 1.d-14
+  INTEGER :: info   !> info code for LAPACK solver
+  INTEGER :: i,j    !> loop variables
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> build a symmetric positive definite matrix in test_A
+  CALL bsm_init(test_A, hbw,n)
+  DO i = 1,n
+    DO j = i,MIN(i+hbw-1,n)
+
+      IF (i.EQ.j) THEN
+        CALL bsm_set(test_A, i,i, ( 2.d0*DBLE(i) ) )  !> i.EQ.j -> 2*i
+      ELSE
+        CALL bsm_set(test_A, i,j, -1.d0)              !> i.NE.j -> -1
+      END IF
+
+    END DO
+  END DO
+
+  !> build test_x = [1..n]
+  CALL vec_init(test_x, n)
+  DO i = 1,n
+    CALL vec_set(test_x, i, DBLE(i))
+  END DO
+
+  !> compute input RHS from test_A*test_x
+  test_b = test_A*test_x
+
+  !> expected solution is test_x
+  expDat = test_x%dat
+
+  !> reset test_x
+  CALL vec_clean(test_x)
+
+  !> compute actual solution, " test_x = test_A^(-1)*test_b "
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, info)
+
+  ! check that solution matches expected result (within tolerance)
+  CALL assertEquals(expDat,test_x%dat, n, toler)
+
+  ! clean up memory allocations
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_VAL
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct solution to linear system of equations (identity)
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix (identity matrix)
+!!  \param test_b   Test right-hand-side vector
+!!  \param test_x   Vector for solution
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expDat   Expected solution
+!!  \param toler    Tolerance on solution (floating point error)
+!!
+!!  This test solves a system of equations, test_A*test_x = test_b, for
+!!  test_x for the case that test_A is the identity matrix. Therefore,
+!!  the expected result is that test_x = test_b.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_IDENT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_IDENT'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  INTEGER, PARAMETER :: hbw=1,n=10
+  DOUBLE PRECISION, DIMENSION(n) :: expDat
+  DOUBLE PRECISION, PARAMETER :: toler = 1.d-14
+  INTEGER :: info   !> info code for LAPACK solver
+  INTEGER :: i,j    !> loop variables
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> build the identity matrix in test_A
+  CALL bsm_init(test_A, hbw,n)
+  DO i = 1,n
+    CALL bsm_set(test_A, i,i, 1.d0)
+  END DO
+
+  !> build test_b = [1..n]
+  CALL vec_init(test_b, n)
+  DO i = 1,n
+    CALL vec_set(test_b, i, DBLE(i))
+  END DO
+
+  !> expected solution is test_b
+  expDat = test_b%dat
+
+  !> compute actual solution, " test_x = test_A^(-1)*test_b "
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, info)
+
+  ! check that solution matches expected result (within tolerance)
+  CALL assertEquals(expDat,test_x%dat, n, toler)
+
+  ! clean up memory allocations
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_IDENT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct solution to linear system of equations (zero)
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix
+!!  \param test_b   Test right-hand-side vector (zeros)
+!!  \param test_x   Vector for solution
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expDat   Expected solution
+!!  \param toler    Tolerance on solution (floating point error)
+!!
+!!  This test solves a system of equations, test_A*test_x = test_b, for
+!!  test_x for the case that test_b is the zero vector. Therefore,
+!!  the expected result is that test_x = test_b = {0}.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_ZERO
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_ZERO'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  INTEGER, PARAMETER :: hbw=3,n=10
+  DOUBLE PRECISION, DIMENSION(n) :: expDat
+  DOUBLE PRECISION, PARAMETER :: toler = 1.d-14
+  INTEGER :: info   !> info code for LAPACK solver
+  INTEGER :: i,j    !> loop variables
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> build a symmetric positive definite matrix in test_A
+  CALL bsm_init(test_A, hbw,n)
+  DO i = 1,n
+    DO j = i,MIN(i+hbw-1,n)
+
+      IF (i.EQ.j) THEN
+        CALL bsm_set(test_A, i,i, ( 2.d0 * DBLE(i) ) )  !> i.EQ.j -> 2*i
+      ELSE
+        CALL bsm_set(test_A, i,j, -1.d0)                !> i.NE.j -> -1
+      END IF
+
+    END DO
+  END DO
+
+  !> initialize test_b (zeros)
+  CALL vec_init(test_b, n)
+
+  !> expected solution is test_b
+  expDat = test_b%dat
+
+  !> compute actual solution, " test_x = test_A^(-1)*test_b "
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, info)
+
+  ! check that solution matches expected result (within tolerance)
+  CALL assertEquals(expDat,test_x%dat, n, toler)
+
+  ! clean up memory allocations
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_ZERO
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for solution in case when decomposition is reused
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param test_A   Test banded symmetric matrix
+!!  \param test_b   Test right-hand-side vector
+!!  \param test_x   Vector for solution
+!!  \param hbw      Half bandwidth of the matrix
+!!  \param n        Number of rows/columns in the matrix
+!!  \param expDat   Expected solution
+!!  \param toler    Tolerance on solution (floating point error)
+!!
+!!  This test solves a system of equations, test_A*test_x = test_b, for
+!!  test_x in the general case for a banded symmetric positive definite
+!!  matrix, test_A. The right-hand-sides, test_b, are computed through
+!!  multiplication of test_A and test_x. That is, the input is
+!!  test_b = test_A*test_x. Therefore, the expected output from the solver
+!!  with inputs test_A and test_b is the original vector test_x. In
+!!  particular, this test solves the system twice. On the second attempt,
+!!  the branch of the linear solver that reuses the decomposition of
+!!  test_A is activated.
+! ------------------------------------------------------------------------
+SUBROUTINE test_linear_solver_DECOMP
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_linear_solver_DECOMP'
+  TYPE(bandSymMatrixT) :: test_A
+  TYPE(vectorT) :: test_b, test_x
+  INTEGER, PARAMETER :: hbw=3,n=10
+  DOUBLE PRECISION, DIMENSION(n) :: expDat
+  DOUBLE PRECISION, PARAMETER :: toler = 1.d-14
+  INTEGER :: info   !> info code for LAPACK solver
+  INTEGER :: i,j    !> loop variables
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> build a symmetric positive definite matrix in test_A
+  CALL bsm_init(test_A, hbw,n)
+  DO i = 1,n
+    DO j = i,MIN(i+hbw-1,n)
+
+      IF (i.EQ.j) THEN
+        CALL bsm_set(test_A, i,i, ( 2.d0 * DBLE(i) ) )  !> i.EQ.j -> 2*i
+      ELSE
+        CALL bsm_set(test_A, i,j, -1.d0)                !> i.NE.j -> -1
+      END IF
+
+    END DO
+  END DO
+
+  !> build test_x = [1..n]
+  CALL vec_init(test_x, n)
+  DO i = 1,n
+    CALL vec_set(test_x, i, DBLE(i))
+  END DO
+
+  !> compute input RHS from test_A*test_x
+  test_b = test_A*test_x
+
+  !> expected solution is test_x
+  expDat = test_x%dat
+
+  !> reset test_x and solve once, " test_x = test_A^(-1)*test_b "
+  CALL vec_clean(test_x)
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, info)
+
+  !> reset test_x and solve again
+  CALL vec_clean(test_x)
+  CALL lin_solve(test_A,test_b, test_x, hbw,n, info)
+
+  ! check that solution matches expected result (within tolerance)
+  CALL assertEquals(expDat,test_x%dat, n, toler)
+
+  ! clean up memory allocations
+  CALL bsm_clean(test_A)
+  CALL vec_clean(test_b)
+  CALL vec_clean(test_x)
+
+END SUBROUTINE test_linear_solver_DECOMP
+
+END MODULE linear_solver_test

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/log_message_control.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Log Message Control module
+! ------------------------------------------------------------------------
+MODULE log_message_control
+USE system_constants    !> Global system constants (for max string length)
+USE log_messages        !> Log/error message and sender codes
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+!> Exported access programs
+PUBLIC :: log_setFileName, log_getFileName, &
+          log_initLogFile, log_closeLogFile, &
+          log_printLogMsg
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+
+INTEGER, PARAMETER :: logFile = 10          !> log file unit
+CHARACTER (LEN=MAXLEN), SAVE :: logName     !> log file name
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for log filename
+!!
+!!  \param fname  String representing filename for problem (without extension)
+!!
+!!  This routine takes the general filename for the problem and appends
+!!  the .log extension to it for the log file.
+! ------------------------------------------------------------------------
+SUBROUTINE log_setFileName (fname)
+  CHARACTER (LEN=*), INTENT(IN) :: fname
+
+  logName = TRIM(fname)//'.log'
+
+END SUBROUTINE log_setFileName
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for log filename
+!!
+!!  \return fname  String representing filename for problem
+! ------------------------------------------------------------------------
+FUNCTION log_getFileName () RESULT(fname)
+  CHARACTER (LEN=MAXLEN) :: fname
+
+  fname = TRIM(logName)
+
+END FUNCTION log_getFileName
+
+
+! ------------------------------------------------------------------------
+!>  \brief Initializer for log file
+!!
+!!  This routine initializes the file for printing log messages. It first
+!!  checks if a log file with the same name as that set in the module
+!!  exists (and deletes it, if present).
+! ------------------------------------------------------------------------
+SUBROUTINE log_initLogFile ()
+  LOGICAL :: fileExists   !> boolean for checking file existence
+
+  !> if file already exists, delete it
+  INQUIRE(FILE=logName, EXIST=fileExists)
+  IF (fileExists) THEN
+    OPEN(UNIT=logFile, FILE=logName, STATUS='OLD')
+    CLOSE(UNIT=logFile, STATUS='DELETE')
+  END IF
+
+  !> open the new file
+  OPEN(UNIT=logFile, FILE=logName)
+
+END SUBROUTINE log_initLogFile
+
+
+! ------------------------------------------------------------------------
+!>  \brief Finalizer for log file
+!!
+!!  This routine first checks if the log file is open, and if so closes
+!!  it.
+! ------------------------------------------------------------------------
+SUBROUTINE log_closeLogFile ()
+  LOGICAL :: fileOpened   !> boolean for checking if file is open
+
+  !> close the log file if it is open
+  INQUIRE(UNIT=logFile, OPENED=fileOpened)
+  IF (fileOpened) THEN
+    CLOSE(UNIT=logFile)
+  END IF
+
+END SUBROUTINE log_closeLogFile
+
+
+! ------------------------------------------------------------------------
+!>  \brief Print log messages to the log file
+!!
+!!  \param msg  Message code (messageT)
+!!  \param sdr  Sender code (senderT)
+!!
+!!  This routine takes the message code and sender code and turns them
+!!  into their associated human readable strings, which it then prints
+!!  to the log file. It makes sure that the log file has been initialized
+!!  and, if not, does not try to print.
+! ------------------------------------------------------------------------
+SUBROUTINE log_printLogMsg (msg, sdr)
+  INTEGER, INTENT(IN) :: msg, sdr
+  LOGICAL :: fileOpened   !> boolean for making sure log file is open
+
+  !> check that log file is open
+  INQUIRE(UNIT=logFile, OPENED=fileOpened)
+
+  !> if so, print the messages
+  IF (fileOpened) THEN
+    WRITE(UNIT=logFile, FMT='(A)') 'Message: '//msg_getMsg(msg)
+    WRITE(UNIT=logFile, FMT='(A)') 'Sender: '//msg_getSdr(sdr)
+    WRITE(UNIT=logFile, FMT=*)
+  END IF
+
+END SUBROUTINE log_printLogMsg
+
+
+END MODULE log_message_control

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/log_message_control_test.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for testing Log Message Control module
+! ------------------------------------------------------------------------
+MODULE log_message_control_test
+USE fruit                 !> Unit testing framework
+USE log_message_control   !> Print log/error messages
+USE log_messages          !> Log/error message and sender codes
+
+IMPLICIT NONE
+
+
+CONTAINS
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for setting of log message filename
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file
+!!  \param expected   Expected filename with extension
+!!
+!!  This test checks that the filename for the log file is parsed and
+!!  appended with the correct extension.
+! ------------------------------------------------------------------------
+SUBROUTINE test_log_setFileName
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_log_setFileName'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  CHARACTER (LEN=*), PARAMETER :: expected = 'testName.log'
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> set the log filename
+  CALL log_setFileName(testName)
+
+  !> check that correct filename is set
+  CALL assertEquals(expected, log_getFileName() )
+
+END SUBROUTINE test_log_setFileName
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for log file initialization
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file
+!!  \param fileExists Boolean for checking existence of log file
+!!
+!!  This test checks that the log file is correctly created.
+! ------------------------------------------------------------------------
+SUBROUTINE test_log_initLogFile
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_log_initLogFile'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  LOGICAL :: fileExists
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> set file name and check if it exists
+  CALL log_setFileName(testName)
+  INQUIRE(FILE=log_getFileName(), EXIST=fileExists)
+
+  !> if file already exists, delete it
+  IF (fileExists) THEN
+    OPEN(UNIT=1, FILE=log_getFileName(), STATUS='OLD')
+    CLOSE(UNIT=1, STATUS='DELETE')
+  END IF
+
+  !> initialize a new file and make sure it exists
+  CALL log_initLogFile()
+  INQUIRE(FILE=log_getFileName(), EXIST=fileExists)
+  CALL assertEquals(.TRUE., fileExists)
+
+END SUBROUTINE test_log_initLogFile
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for log file finalization
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file
+!!  \param fileOpened Boolean for checking whether file is open
+!!
+!!  This test checks that the log file is correctly closed.
+! ------------------------------------------------------------------------
+SUBROUTINE test_log_closeLogFile
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_log_closeLogFile'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  LOGICAL :: fileOpened
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> set file name and initialize
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> close log file and make sure file is not associated with a unit
+  CALL log_closeLogFile()
+  INQUIRE(FILE=log_getFileName(), OPENED=fileOpened)
+  CALL assertEquals(.FALSE., fileOpened)
+
+END SUBROUTINE test_log_closeLogFile
+
+
+! ------------------------------------------------------------------------
+!>  \test Test of log message printing
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file
+!!  \param fileOpened Boolean for checking whether file is open
+!!
+!!  This test checks for correct printing of a log message to the log
+!!  file.
+! ------------------------------------------------------------------------
+SUBROUTINE test_log_printLogMsg
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_log_printLogMsg'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  CHARACTER (LEN=*), PARAMETER :: expMsg = 'Message: Failed to allocate memory'
+  CHARACTER (LEN=*), PARAMETER :: expSdr = 'Sender: Body Force Reader'
+  CHARACTER (LEN=MAXLEN) :: actMsg, actSdr
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> set file name and initialize
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> write log message
+  CALL log_printLogMsg(ALLOC, BFCRDR)
+
+  !> close log file
+  CALL log_closeLogFile()
+
+  !> open log file and read for test
+  OPEN(UNIT=1, FILE=log_getFileName())
+  READ(UNIT=1, FMT='(A)', ERR=1, END=1) actMsg
+  READ(UNIT=1, FMT='(A)', ERR=1, END=1) actSdr
+1 CLOSE(UNIT=1)   !> label must be used here to avoid file reading errors
+
+  !> check that correct messages were written
+  CALL assertEquals(expMsg, actMsg)
+  CALL assertEquals(expSdr, actSdr)
+
+END SUBROUTINE test_log_printLogMsg
+
+END MODULE log_message_control_test

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/log_messages.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module defining Log Messages
+! ------------------------------------------------------------------------
+MODULE log_messages
+USE system_constants    !> Global system constants (for max string length)
+
+IMPLICIT NONE
+
+
+
+! ************************************************************************
+! EXPORTED DATA TYPES
+! ************************************************************************
+
+
+!> \brief Enumerated type messageT
+ENUM, BIND(C)
+  ENUMERATOR :: OK=1, ALLOC=2, DIMEN=3, EXCEED=4, EXISTS=5, FORMT=6
+  ENUMERATOR :: POSIT=7, POSDEF=8, SZE=9, TYP=10
+END ENUM ! messageT
+
+!> \brief Enumerated type senderT
+ENUM, BIND(C)
+  ENUMERATOR :: BFCRDR=1, BNDDAT=2, BNDRDR=3, BSYMAT=4, CNSMAT=5, DMNRDR=6
+  ENUMERATOR :: DNSMAT=7, FLDDAT=8, ICTRDR=9, ICVRDR=10, KBCRDR=11
+  ENUMERATOR :: LINSLV=12, MTLDAT=13, MTLRDR=14, NBCRDR=15, TNSWTR=16
+  ENUMERATOR :: VECTOR=17, VECWTR=18
+END ENUM ! senderT
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! EXPORTED ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Function for getting the string message corresponding to a messageT code
+!!
+!!  \param  code Message code (messageT)
+!!  \return exc  Human readable string representing exception
+! ------------------------------------------------------------------------
+FUNCTION msg_getMsg (code) RESULT(exc)
+  INTEGER, INTENT(IN) :: code
+  CHARACTER(LEN=MAXLEN) :: exc
+
+  SELECT CASE (code)
+    CASE (ALLOC)
+      exc = 'Failed to allocate memory'
+    CASE (DIMEN)
+      exc = 'Size of data structure does not match size of target data structure'
+    CASE (EXCEED)
+      exc = 'Data value exceeds a defined minimum or maximum'
+    CASE (EXISTS)
+      exc = 'File does not exist'
+    CASE (FORMT)
+      exc = 'File data is not in expected format'
+    CASE (POSDEF)
+      exc = 'Matrix is not positive definite'
+    CASE (POSIT)
+      exc = 'Index exceeds size of data structure'
+    CASE (SZE)
+      exc = 'Specified size of data structure exceeds minimum or maximum allowable size'
+    CASE (TYP)
+      exc = 'Specified material property does not correspond to material type'
+    CASE DEFAULT
+      exc = 'Could not find log message code'
+  END SELECT
+
+END FUNCTION msg_getMsg
+
+
+! ------------------------------------------------------------------------
+!>  \brief Function for getting the string message corresponding to a senderT code
+!!
+!!  \param  code Message code (senderT)
+!!  \return sdr  Human readable string representing sender
+! ------------------------------------------------------------------------
+FUNCTION msg_getSdr (code) RESULT(sdr)
+  INTEGER, INTENT(IN) :: code
+  CHARACTER(LEN=MAXLEN) :: sdr
+
+  SELECT CASE (code)
+    CASE (BFCRDR)
+      sdr = 'Body Force Reader'
+    CASE (BNDDAT)
+      sdr = 'Boundary Data'
+    CASE (BNDRDR)
+      sdr = 'Boundary File Reader'
+    CASE (BSYMAT)
+      sdr = 'Banded Symmetric Matrix'
+    CASE (CNSMAT)
+      sdr = 'Constitutive Matrix'
+    CASE (DMNRDR)
+      sdr = 'Domain File Reader'
+    CASE (DNSMAT)
+      sdr = 'Dense Matrix'
+    CASE (FLDDAT)
+      sdr = 'Field Data'
+    CASE (ICTRDR)
+      sdr = 'Initial Tensor Field Reader'
+    CASE (ICVRDR)
+      sdr = 'Initial Vector Field Reader'
+    CASE (KBCRDR)
+      sdr = 'Kinematic BC Reader'
+    CASE (LINSLV)
+      sdr = 'Linear Solver'
+    CASE (MTLDAT)
+      sdr = 'Material Property Data'
+    CASE (MTLRDR)
+      sdr = 'Material File Reader'
+    CASE (NBCRDR)
+      sdr = 'Natural BC Reader'
+    CASE (TNSWTR)
+      sdr = 'Tensor Field Writer'
+    CASE (VECTOR)
+      sdr = 'Vector Data Type'
+    CASE (VECWTR)
+      sdr = 'Vector Field Writer'
+    CASE DEFAULT
+      sdr = 'Could not find sender code'
+  END SELECT
+
+END FUNCTION msg_getSdr
+
+END MODULE log_messages

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/log_messages_test.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for testing Log Messages module
+! ------------------------------------------------------------------------
+MODULE log_messages_test
+USE fruit             !> Unit testing framework
+USE system_constants  !> Global system constants (for max string length)
+USE log_messages      !> Log/error message and sender codes
+
+IMPLICIT NONE
+
+
+CONTAINS
+
+
+! ************************************************************************
+! TEST MESSAGE CODES
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for ALLOC message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_ALLOC_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_ALLOC_message_code'
+  INTEGER, PARAMETER :: code = ALLOC
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Failed to allocate memory'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_ALLOC_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for DIMEN message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_DIMEN_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_DIMEN_message_code'
+  INTEGER, PARAMETER :: code = DIMEN
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Size of data structure does not match size of target data structure'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_DIMEN_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for EXCEED message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_EXCEED_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_EXCEED_message_code'
+  INTEGER, PARAMETER :: code = EXCEED
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Data value exceeds a defined minimum or maximum'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_EXCEED_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for EXISTS message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_EXISTS_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_EXISTS_message_code'
+  INTEGER, PARAMETER :: code = EXISTS
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'File does not exist'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_EXISTS_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for FORMT message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_FORMT_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_FORMT_message_code'
+  INTEGER, PARAMETER :: code = FORMT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'File data is not in expected format'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_FORMT_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_POSIT_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_POSIT_message_code'
+  INTEGER, PARAMETER :: code = POSIT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Index exceeds size of data structure'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_POSIT_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSDEF message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_POSDEF_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_POSDEF_message_code'
+  INTEGER, PARAMETER :: code = POSDEF
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Matrix is not positive definite'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_POSDEF_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for SZE message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_SZE_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_SZE_message_code'
+  INTEGER, PARAMETER :: code = SZE
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Specified size of data structure exceeds minimum or maximum allowable size'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_SZE_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for TYP message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_TYP_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_TYP_message_code'
+  INTEGER, PARAMETER :: code = TYP
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Specified material property does not correspond to material type'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_TYP_message_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for unexpected message code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Message code (not in enumerated list of message codes)
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for a
+!!  non-existent message code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_unexpected_message_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_unexpected_message_code'
+  INTEGER, PARAMETER :: code = 100
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Could not find log message code'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getMsg(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_unexpected_message_code
+
+
+
+! ************************************************************************
+! TEST SENDER CODES
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for BFCRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_BFCRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_BFCRDR_sender_code'
+  INTEGER, PARAMETER :: code = BFCRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Body Force Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_BFCRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for BNDDAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_BNDDAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_BNDDAT_sender_code'
+  INTEGER, PARAMETER :: code = BNDDAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Boundary Data'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_BNDDAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for BNDRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_BNDRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_BNDRDR_sender_code'
+  INTEGER, PARAMETER :: code = BNDRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Boundary File Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_BNDRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for BSYMAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_BSYMAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_BSYMAT_sender_code'
+  INTEGER, PARAMETER :: code = BSYMAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Banded Symmetric Matrix'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_BSYMAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for CNSMAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_CNSMAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_CNSMAT_sender_code'
+  INTEGER, PARAMETER :: code = CNSMAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Constitutive Matrix'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_CNSMAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for DMNRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_DMNRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_DMNRDR_sender_code'
+  INTEGER, PARAMETER :: code = DMNRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Domain File Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_DMNRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for DNSMAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_DNSMAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_DNSMAT_sender_code'
+  INTEGER, PARAMETER :: code = DNSMAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Dense Matrix'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_DNSMAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for FLDDAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_FLDDAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_FLDDAT_sender_code'
+  INTEGER, PARAMETER :: code = FLDDAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Field Data'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_FLDDAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for ICTRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_ICTRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_ICTRDR_sender_code'
+  INTEGER, PARAMETER :: code = ICTRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Initial Tensor Field Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_ICTRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for ICVRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_ICVRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_ICVRDR_sender_code'
+  INTEGER, PARAMETER :: code = ICVRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Initial Vector Field Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_ICVRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for KBCRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_KBCRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_KBCRDR_sender_code'
+  INTEGER, PARAMETER :: code = KBCRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Kinematic BC Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_KBCRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for LINSLV sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_LINSLV_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_LINSLV_sender_code'
+  INTEGER, PARAMETER :: code = LINSLV
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Linear Solver'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_LINSLV_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for MTLDAT sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_MTLDAT_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_MTLDAT_sender_code'
+  INTEGER, PARAMETER :: code = MTLDAT
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Material Property Data'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_MTLDAT_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for MTLRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_MTLRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_MTLRDR_sender_code'
+  INTEGER, PARAMETER :: code = MTLRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Material File Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_MTLRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for NBCRDR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_NBCRDR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_NBCRDR_sender_code'
+  INTEGER, PARAMETER :: code = NBCRDR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Natural BC Reader'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_NBCRDR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for TNSWTR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_TNSWTR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_TNSWTR_sender_code'
+  INTEGER, PARAMETER :: code = TNSWTR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Tensor Field Writer'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_TNSWTR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for VECTOR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_VECTOR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_VECTOR_sender_code'
+  INTEGER, PARAMETER :: code = VECTOR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Vector Data Type'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_VECTOR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for VECWTR sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for the
+!!  corresponding sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_VECWTR_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_VECWTR_sender_code'
+  INTEGER, PARAMETER :: code = VECWTR
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Vector Field Writer'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_VECWTR_sender_code
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for unexpected sender code
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param code       Sender code (not in enumerated list of sender codes)
+!!  \param expected   Expected message
+!!  \param actual     Actual message
+!!
+!!  This test checks that the correct message is returned for a
+!!  non-existent sender code.
+! ------------------------------------------------------------------------
+SUBROUTINE test_unexpected_sender_code
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_unexpected_sender_code'
+  INTEGER, PARAMETER :: code = 100
+  CHARACTER(LEN=MAXLEN), PARAMETER :: expected = 'Could not find sender code'
+  CHARACTER(LEN=MAXLEN) :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  actual = msg_getSdr(code)
+
+  CALL assertEquals(TRIM(expected),TRIM(actual))
+
+END SUBROUTINE test_unexpected_sender_code
+
+
+END MODULE log_messages_test

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/material_data.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for Material Property Data
+! ------------------------------------------------------------------------
+MODULE material_data
+USE system_constants        !> Global system constants (for min/max prop values)
+USE log_message_control     !> Print log/error messages
+USE log_messages            !> Log/error codes and messages
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: mtl_init, mtl_clean, &
+          mtl_numMtl, &
+          mtl_getEmod, mtl_setEmod, &
+          mtl_getPois, mtl_setPois, &
+          mtl_getDens, mtl_setDens
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+!> sender code for this module
+INTEGER, PARAMETER :: sdr = MTLDAT
+
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+
+!> enumerated type indicating material model
+ENUM, BIND(C)
+  ENUMERATOR :: linear_elastic
+END ENUM
+
+
+! ------------------------------------------------------------------------
+!>  \brief Material Type
+!!
+!!  \param num    Hash value indicating material number
+!!  \param typ    Material model type (from enumerated type in this module)
+!!  \param emod   Elastic modulus
+!!  \param nu     Poisson's ratio
+!!  \param rho    Density
+! ------------------------------------------------------------------------
+TYPE materialT
+  INTEGER :: num
+  INTEGER :: typ
+  DOUBLE PRECISION :: emod, nu, rho
+END TYPE materialT
+
+
+
+! ************************************************************************
+! STATE VARIABLES
+! ************************************************************************
+
+TYPE(materialT), ALLOCATABLE :: materials(:)
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to constructor for module state variables
+INTERFACE mtl_init
+  MODULE PROCEDURE mtl_init_
+  MODULE PROCEDURE mtl_init_exc_
+END INTERFACE mtl_init
+
+!> \brief Interface to destructor for module state variables
+INTERFACE mtl_clean
+  MODULE PROCEDURE mtl_clean_
+END INTERFACE mtl_clean
+
+!> \brief Interface for number of materials
+INTERFACE mtl_numMtl
+  MODULE PROCEDURE mtl_num_mtl_
+END INTERFACE mtl_numMtl
+
+!> \brief Interface for getter for elastic modulus
+INTERFACE mtl_getEmod
+  MODULE PROCEDURE mtl_get_emod_
+  MODULE PROCEDURE mtl_get_emod_exc_
+END INTERFACE mtl_getEmod
+
+!> \brief Interface for setter for elastic modulus
+INTERFACE mtl_setEmod
+  MODULE PROCEDURE mtl_set_emod_
+  MODULE PROCEDURE mtl_set_emod_exc_
+END INTERFACE mtl_setEmod
+
+!> \brief Interface for getter for Poisson's ratio
+INTERFACE mtl_getPois
+  MODULE PROCEDURE mtl_get_pois_
+  MODULE PROCEDURE mtl_get_pois_exc_
+END INTERFACE mtl_getPois
+
+!> \brief Interface for setter for Poisson's ratio
+INTERFACE mtl_setPois
+  MODULE PROCEDURE mtl_set_pois_
+  MODULE PROCEDURE mtl_set_pois_exc_
+END INTERFACE mtl_setPois
+
+!> \brief Interface for getter for density
+INTERFACE mtl_getDens
+  MODULE PROCEDURE mtl_get_dens_
+  MODULE PROCEDURE mtl_get_dens_exc_
+END INTERFACE mtl_getDens
+
+!> \brief Interface for setter for density
+INTERFACE mtl_setDens
+  MODULE PROCEDURE mtl_set_dens_
+  MODULE PROCEDURE mtl_set_dens_exc_
+END INTERFACE mtl_setDens
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Constructor for Material Data state variables (non-exception checking)
+!!
+!!  \param  nmtl  Number of materials
+!!
+!!  This routine allocates memory and initializes the state variable that
+!!  contains the set of material property information for the problem.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_init_ (nmtl)
+  INTEGER, INTENT(IN) :: nmtl
+  INTEGER :: imtl     !> loop variable
+
+  !> only reallocate if new dimension does not match existing dimension
+  IF (mtl_numMtl().NE.nmtl) THEN
+
+    !> ensure state variable is clear
+    CALL mtl_clean()
+
+    !> allocate memory for material data state variable
+    ALLOCATE(materials(nmtl))
+
+  END IF
+
+  !> initialize state variable
+  DO imtl = 1,nmtl
+    materials(imtl)%num   = imtl
+    materials(imtl)%typ   = linear_elastic
+    materials(imtl)%emod  = 0.d0
+    materials(imtl)%nu    = 0.d0
+    materials(imtl)%rho   = 0.d0
+  END DO
+
+END SUBROUTINE mtl_init_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Constructor for Material Data state variables (exception checking)
+!!
+!!  \param      nmtl  Number of materials
+!!  \param      exc   Error code
+!!
+!!  \exception  ALLOC Memory allocation for state variable failed
+!!  \exception  SZE   Specified number of materials is invalid
+!!
+!!  This routine allocates memory and initializes the state variable that
+!!  contains the set of material property information for the problem.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_init_exc_ (nmtl, exc)
+  INTEGER, INTENT(IN) :: nmtl
+  INTEGER, INTENT(OUT) :: exc
+  INTEGER :: e        !> error code for allocation
+  INTEGER :: imtl     !> loop variable
+
+  !> ensure that specified number of materials is valid
+  IF (nmtl.LT.1 .OR. nmtl.GT.MAX_MATERIALS) THEN
+    exc=SZE
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+  ELSE
+    exc=OK
+  END IF
+
+  !> only reallocate if new dimension does not match existing dimension
+  IF (mtl_numMtl().NE.nmtl) THEN
+
+    !> ensure state variable is clear
+    CALL mtl_clean()
+
+    !> allocate memory for material data state variable
+    ALLOCATE(materials(nmtl), STAT=e)
+
+    !> ensure that memory allocation was successful
+    IF (e.NE.0) THEN
+      exc=ALLOC
+      CALL log_printLogMsg(exc,sdr)
+      RETURN
+    ELSE
+      exc=OK
+    END IF
+
+  END IF
+
+  !> initialize state variable
+  DO imtl = 1,nmtl
+    materials(imtl)%num   = imtl
+    materials(imtl)%typ   = linear_elastic
+    materials(imtl)%emod  = 0.d0
+    materials(imtl)%nu    = 0.d0
+    materials(imtl)%rho   = 0.d0
+  END DO
+
+END SUBROUTINE mtl_init_exc_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Destructor for Material Data state variables
+!!
+!!  This routine clears the memory allocated to the state variable that
+!!  contains the set of material property information for the problem.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_clean_ ()
+
+  IF (ALLOCATED(materials)) DEALLOCATE(materials)
+
+END SUBROUTINE mtl_clean_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for number of materials
+!!
+!!  \return nmtl  Number of materials
+!!
+!!  This routine determines the number of data entries that have been
+!!  allocated for material data. It does not check that the material data
+!!  has been populated (i.e. changed from initial zero values).
+! ------------------------------------------------------------------------
+FUNCTION mtl_num_mtl_ () RESULT(nmtl)
+  INTEGER :: nmtl
+
+  !> if data is initialized, return number of materials
+  IF (ALLOCATED(materials)) THEN
+    nmtl = SIZE(materials)
+  ELSE
+    nmtl = 0   !> if not initialized, there are no materials
+  END IF
+
+END FUNCTION mtl_num_mtl_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for elastic modulus (non-exception checking)
+!!
+!!  \param  i     Material number
+!!  \return emod  Elastic modulus
+!!
+!!  This routine determines the value of the elastic modulus for material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_emod_ (i) RESULT(emod)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION :: emod
+
+  emod = materials(i)%emod
+
+END FUNCTION mtl_get_emod_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for elastic modulus (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  exc   Error code
+!!
+!!  \return emod  Elastic modulus
+!!
+!!  \exception  POSIT The material number is not in [1..mtl_numMtls()]
+!!  \exception  TYP   The material type is not 'linear_elastic'
+!!
+!!  This routine determines the value of the elastic modulus for material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_emod_exc_ (i, exc) RESULT(emod)
+  INTEGER, INTENT(IN) :: i
+  INTEGER, INTENT(OUT) :: exc
+  DOUBLE PRECISION :: emod
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    emod = 0.d0
+    RETURN
+
+  !> check that the material type is linear elastic
+  !! (otherwise, elastic modulus may not be appropriate)
+  ELSE IF ( materials(i)%typ .NE. linear_elastic ) THEN
+    exc=TYP
+    CALL log_printLogMsg(exc,sdr)
+    emod = 0.d0
+    RETURN
+
+  ELSE
+    exc=OK
+  END IF
+
+  !> call non-exception version
+  emod = mtl_getEmod(i)
+
+END FUNCTION mtl_get_emod_exc_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for elastic modulus (non-exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  emod  Elastic modulus
+!!
+!!  This routine sets the value of the elastic modulus for material
+!!  number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_emod_ (i,emod)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: emod
+
+  materials(i)%emod = emod
+
+END SUBROUTINE mtl_set_emod_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for elastic modulus (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  emod  Elastic modulus
+!!  \param  exc   Error code
+!!
+!!  \exception  EXCEED  The value of elastic modulus exceeds defined limits
+!!  \exception  POSIT   The material number is not in [1..mtl_numMtls()]
+!!  \exception  TYP     The material type is not 'linear_elastic'
+!!
+!!  This routine sets the value of the elastic modulus for material
+!!  number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_emod_exc_ (i,emod, exc)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: emod
+  INTEGER, INTENT(OUT) :: exc
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  !> check that the material type is linear elastic
+  !! (otherwise, elastic modulus may not be appropriate)
+  ELSE IF ( materials(i)%typ .NE. linear_elastic ) THEN
+    exc=TYP
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  !> check that the value of elastic modulus is within prescribed limits
+  ELSE IF (emod.LT.E_MIN .OR. emod.GT.E_MAX) THEN
+    exc=EXCEED
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  ELSE
+    exc=OK
+
+  END IF
+
+  !> call non-exception version
+  CALL mtl_setEmod(i,emod)
+
+END SUBROUTINE mtl_set_emod_exc_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for Poisson's ratio (non-exception checking)
+!!
+!!  \param  i   Material number
+!!  \return nu  Poisson's ratio
+!!
+!!  This routine determines the value of Poisson's ratio for material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_pois_ (i) RESULT(nu)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION :: nu
+
+  nu = materials(i)%nu
+
+END FUNCTION mtl_get_pois_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for Poisson's ratio (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  exc   Error code
+!!
+!!  \return nu    Poisson's ratio
+!!
+!!  \exception  POSIT The material number is not in [1..mtl_numMtls()]
+!!  \exception  TYP   The material type is not 'linear_elastic'
+!!
+!!  This routine determines the value of Poisson's ratio for material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_pois_exc_ (i, exc) RESULT(nu)
+  INTEGER, INTENT(IN) :: i
+  INTEGER, INTENT(OUT) :: exc
+  DOUBLE PRECISION :: nu
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    nu = 0.d0
+    RETURN
+
+  !> check that the material type is linear elastic
+  !! (otherwise, Poisson's ratio may not be appropriate)
+  ELSE IF ( materials(i)%typ .NE. linear_elastic ) THEN
+    exc=TYP
+    CALL log_printLogMsg(exc,sdr)
+    nu = 0.d0
+    RETURN
+
+  ELSE
+    exc=OK
+  END IF
+
+  !> call non-exception version
+  nu = mtl_getPois(i)
+
+END FUNCTION mtl_get_pois_exc_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for Poisson's ratio (non-exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  nu    Poisson's ratio
+!!
+!!  This routine sets the value of Poisson's ratio for material
+!!  number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_pois_ (i,nu)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: nu
+
+  materials(i)%nu = nu
+
+END SUBROUTINE mtl_set_pois_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for Poisson's ratio (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  nu    Poisson's ratio
+!!  \param  exc   Error code
+!!
+!!  \exception  EXCEED  The value of Poisson's ratio exceeds defined limits
+!!  \exception  POSIT   The material number is not in [1..mtl_numMtls()]
+!!  \exception  TYP     The material type is not 'linear_elastic'
+!!
+!!  This routine sets the value of Poisson's ratio for material
+!!  number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_pois_exc_ (i,nu, exc)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: nu
+  INTEGER, INTENT(OUT) :: exc
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  !> check that the material type is linear elastic
+  !! (otherwise, elastic modulus may not be appropriate)
+  ELSE IF ( materials(i)%typ .NE. linear_elastic ) THEN
+    exc=TYP
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  !> check that the value of elastic modulus is within prescribed limits
+  ELSE IF (nu.LT.NU_MIN .OR. nu.GT.NU_MAX) THEN
+    exc=EXCEED
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  ELSE
+    exc=OK
+
+  END IF
+
+  !> call non-exception version
+  CALL mtl_setPois(i,nu)
+
+END SUBROUTINE mtl_set_pois_exc_
+
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for density (non-exception checking)
+!!
+!!  \param  i     Material number
+!!  \return rho   Density
+!!
+!!  This routine determines the value of the density of material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_dens_ (i) RESULT(rho)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION :: rho
+
+  rho = materials(i)%rho
+
+END FUNCTION mtl_get_dens_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Getter for density (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  exc   Error code
+!!
+!!  \return rho    Density
+!!
+!!  \exception  POSIT The material number is not in [1..mtl_numMtls()]
+!!
+!!  This routine determines the value of the density of material
+!!  number i.
+! ------------------------------------------------------------------------
+FUNCTION mtl_get_dens_exc_ (i, exc) RESULT(rho)
+  INTEGER, INTENT(IN) :: i
+  INTEGER, INTENT(OUT) :: exc
+  DOUBLE PRECISION :: rho
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    rho = 0.d0
+    RETURN
+  ELSE
+    exc=OK
+  END IF
+
+  !> call non-exception version
+  rho = mtl_getDens(i)
+
+END FUNCTION mtl_get_dens_exc_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for density (non-exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  rho   Density
+!!
+!!  This routine sets the value of the density of material number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_dens_ (i,rho)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: rho
+
+  materials(i)%rho = rho
+
+END SUBROUTINE mtl_set_dens_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Setter for density (exception checking)
+!!
+!!  \param  i     Material number
+!!  \param  rho   Density
+!!  \param  exc   Error code
+!!
+!!  \exception  EXCEED  The value of the density exceeds defined limits
+!!  \exception  POSIT   The material number is not in [1..mtl_numMtls()]
+!!
+!!  This routine sets the value of the density of material number i.
+! ------------------------------------------------------------------------
+SUBROUTINE mtl_set_dens_exc_ (i,rho, exc)
+  INTEGER, INTENT(IN) :: i
+  DOUBLE PRECISION, INTENT(IN) :: rho
+  INTEGER, INTENT(OUT) :: exc
+
+  !> check that the index is within the bounds of the material list
+  IF ( i.LT.1 .OR. i.GT.mtl_numMtl() ) THEN
+    exc=POSIT
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  !> check that the value of density is within prescribed limits
+  ELSE IF (rho.LT.RHO_MIN .OR. rho.GT.RHO_MAX) THEN
+    exc=EXCEED
+    CALL log_printLogMsg(exc,sdr)
+    RETURN
+
+  ELSE
+    exc=OK
+
+  END IF
+
+  !> call non-exception version
+  CALL mtl_setDens(i,rho)
+
+END SUBROUTINE mtl_set_dens_exc_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE material_data

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/material_data_test.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for testing Material Property Data module
+! ------------------------------------------------------------------------
+MODULE material_data_test
+USE fruit                   !> Unit testing framework
+USE system_constants        !> Global constants
+USE log_message_control     !> Printing log/error messages
+USE log_messages            !> Log/error codes
+USE material_data           !> Material Property Data module
+
+IMPLICIT NONE
+
+
+CONTAINS
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for OK exception message on allocation
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!
+!!  This test initializes the material_data module and makes sure that the
+!!  exception message is OK (i.e. allocation did not fail).
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_allocation_MSG
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_allocation_MSG'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg=OK
+  INTEGER :: actMsg
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log message file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize module and check the exception
+  CALL mtl_init(nmtl, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize log file and deallocate module
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_allocation_MSG
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for SZE exception message on allocation
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!
+!!  This test attempts to initialize the material data module with invalid
+!!  size parameters and verifies that the correct exception is returned.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_allocation_SZE
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_allocation_SZE'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: expMsg=SZE
+  INTEGER :: actMsg
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> try to initialize with nmtl=0
+  CALL mtl_init(0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> ensure module is reset
+  CALL mtl_clean()
+
+  !> try to initialize with nmtl=MAX_MATERIALS+1
+  CALL mtl_init(MAX_MATERIALS+1, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize log file and deallocate module
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_allocation_SZE
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for number of materials when module is not initialized
+!!
+!!  \param unit_name    Name of unit test (for FRUIT)
+!!  \param expected     Expected number of materials
+!!  \param actual       Actual number of materials
+!!
+!!  This test makes sure that the number of materials is returned as 0
+!!  when the module is not initialized
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_num_mtl_not_allocated
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_num_mtl_not_allocated'
+  INTEGER, PARAMETER :: expected = 0
+  INTEGER :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> check number of materials
+  actual = mtl_numMtl()
+  CALL assertEquals(expected, actual)
+
+END SUBROUTINE test_mtl_num_mtl_not_allocated
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for number of materials when module is initialized
+!!
+!!  \param unit_name    Name of unit test (for FRUIT)
+!!  \param expected     Expected number of materials
+!!  \param actual       Actual number of materials
+!!
+!!  This test makes sure that the correct number of materials is returned.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_num_mtl_allocated
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_num_mtl_allocated'
+  INTEGER, PARAMETER :: expected = 10
+  INTEGER :: actual
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize the module
+  CALL mtl_init(expected)
+
+  !> check number of materials
+  actual = mtl_numMtl()
+  CALL assertEquals(expected, actual)
+
+  !> deallocate the module
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_num_mtl_allocated
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_getEmod
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param emod       Dummy variable for get function return
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_emod_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_emod_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION :: emod
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to get beyond last material
+  emod = mtl_getEmod(nmtl+1, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to get before first material
+  emod = mtl_getEmod(0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_emod_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_setEmod
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param emod       Dummy variable for set routine input
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_emod_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_emod_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION, PARAMETER :: emod = 3.d0
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set beyond last material
+  CALL mtl_setEmod(nmtl+1,emod, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set before first material
+  CALL mtl_setEmod(0,emod, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_emod_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for EXCEED exception from mtl_setEmod
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!
+!!  This test checks that an EXCEED exception is returned when the
+!!  specified input is not within the range defined in the System
+!!  Constants module.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_emod_EXCEED
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_emod_EXCEED'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  INTEGER, PARAMETER :: expMsg = EXCEED
+  INTEGER :: actMsg
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set below min value
+  CALL mtl_setEmod(i,E_MIN-1.d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set above max value
+  CALL mtl_setEmod(i,E_MAX+1.d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_emod_EXCEED
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct value getting and setting in mtl_getEmod and mtl_setEmod
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expVal     Expected data value
+!!  \param actVal     Actual data value
+!!
+!!  This test checks that the correct value is set using the mtl_setEmod
+!!  access program and returned from the mtl_getEmod access program.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_set_emod_VAL
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_set_emod_VAL'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  DOUBLE PRECISION, PARAMETER :: expVal = 3.d0
+  DOUBLE PRECISION :: actVal
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> set the value of the material
+  CALL mtl_setEmod(i,expVal)
+
+  !> get the value using the access program
+  actVal = mtl_getEmod(i)
+  CALL assertEquals(expVal,actVal)
+
+  !> deallocate the materials
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_set_emod_VAL
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_getPois
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param nu         Dummy variable for get function return
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_pois_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_pois_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION :: nu
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to get beyond last material
+  nu = mtl_getPois(nmtl+1, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to get before first material
+  nu = mtl_getPois(0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_pois_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_setPois
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param nu         Dummy variable for set routine input
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_pois_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_pois_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION, PARAMETER :: nu = 0.25d0
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set beyond last material
+  CALL mtl_setPois(nmtl+1,nu, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set before first material
+  CALL mtl_setPois(0,nu, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_pois_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for EXCEED exception from mtl_setPois
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!
+!!  This test checks that an EXCEED exception is returned when the
+!!  specified input is not within the range defined in the System
+!!  Constants module.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_pois_EXCEED
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_pois_EXCEED'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  INTEGER, PARAMETER :: expMsg = EXCEED
+  INTEGER :: actMsg
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set below min value
+  CALL mtl_setPois(i,NU_MIN-0.1d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set above max value
+  CALL mtl_setPois(i,NU_MAX+0.1d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_pois_EXCEED
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct value getting and setting in mtl_getPois and mtl_setPois
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expVal     Expected data value
+!!  \param actVal     Actual data value
+!!
+!!  This test checks that the correct value is set using the mtl_setPois
+!!  access program and returned from the mtl_getPois access program.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_set_pois_VAL
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_set_pois_VAL'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  DOUBLE PRECISION, PARAMETER :: expVal = 0.25d0
+  DOUBLE PRECISION :: actVal
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> set the value of the material
+  CALL mtl_setPois(i,expVal)
+
+  !> get the value using the access program
+  actVal = mtl_getPois(i)
+  CALL assertEquals(expVal,actVal)
+
+  !> deallocate the materials
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_set_pois_VAL
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_getDens
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param rho        Dummy variable for get function return
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_dens_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_dens_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION :: rho
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to get beyond last material
+  rho = mtl_getDens(nmtl+1, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to get before first material
+  rho = mtl_getDens(0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_dens_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for POSIT exception from mtl_setDens
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!  \param rho        Dummy variable for set routine input
+!!
+!!  This test checks that a POSIT exception is returned when the requested
+!!  location is not inside the material list.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_dens_POSIT
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_dens_POSIT'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: expMsg = POSIT
+  INTEGER :: actMsg
+  DOUBLE PRECISION, PARAMETER :: rho = 2.d3
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set beyond last material
+  CALL mtl_setDens(nmtl+1,rho, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set before first material
+  CALL mtl_setDens(0,rho, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_dens_POSIT
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for EXCEED exception from mtl_setDens
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param testName   Filename for log file (required for exceptions)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expMsg     Expected log message
+!!  \param actMsg     Actual log message
+!!
+!!  This test checks that an EXCEED exception is returned when the
+!!  specified input is not within the range defined in the System
+!!  Constants module.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_set_dens_EXCEED
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_set_dens_EXCEED'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'testName'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  INTEGER, PARAMETER :: expMsg = EXCEED
+  INTEGER :: actMsg
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize log file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> try to set below min value
+  CALL mtl_setDens(i,RHO_MIN-1.d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> try to set above max value
+  CALL mtl_setDens(i,RHO_MAX+1.d0, actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> finalize the log file and deallocate the material list
+  CALL log_closeLogFile()
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_set_dens_EXCEED
+
+
+! ------------------------------------------------------------------------
+!>  \test Test for correct value getting and setting in mtl_getDens and mtl_setDens
+!!
+!!  \param unit_name  Name of unit test (for FRUIT)
+!!  \param nmtl       Number of materials
+!!  \param i          Material number for test location
+!!  \param expVal     Expected data value
+!!  \param actVal     Actual data value
+!!
+!!  This test checks that the correct value is set using the mtl_setDens
+!!  access program and returned from the mtl_getDens access program.
+! ------------------------------------------------------------------------
+SUBROUTINE test_mtl_get_set_dens_VAL
+  CHARACTER (LEN=*), PARAMETER :: unit_name = 'test_mtl_get_set_dens_VAL'
+  INTEGER, PARAMETER :: nmtl=10
+  INTEGER, PARAMETER :: i=3
+  DOUBLE PRECISION, PARAMETER :: expVal = 2.d3
+  DOUBLE PRECISION :: actVal
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name)
+
+  !> initialize materials
+  CALL mtl_init(nmtl)
+
+  !> set the value of the material
+  CALL mtl_setDens(i,expVal)
+
+  !> get the value using the access program
+  actVal = mtl_getDens(i)
+  CALL assertEquals(expVal,actVal)
+
+  !> deallocate material data
+  CALL mtl_clean()
+
+END SUBROUTINE test_mtl_get_set_dens_VAL
+
+
+END MODULE material_data_test

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/material_model.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for Material Models
+! ------------------------------------------------------------------------
+MODULE material_model
+USE dense_matrix_def    !> Dense Matrix data type
+USE vector_def          !> Vector data type
+USE constitutive        !> Constitutive Matrix module
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported interfaces
+PUBLIC :: linearElastic
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+
+!> \brief Interface to linear elastic material model
+INTERFACE linearElastic
+  MODULE PROCEDURE linear_elastic_2d_plane_strain_
+END INTERFACE linearElastic
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Linear elastic material model (2-d, plane strain)
+!!
+!!  \param  emod  Elastic modulus
+!!  \param  nu    Poisson's ratio
+!!  \param  deps  Strain increment
+!!  \param  dsig  Stress increment
+!!
+!!  This routine computes the incremental stress given the incremental
+!!  strain for a linear elastic material.
+! ------------------------------------------------------------------------
+SUBROUTINE linear_elastic_2d_plane_strain_ (emod,nu, deps, dsig)
+  DOUBLE PRECISION, INTENT(IN) :: emod,nu
+  TYPE(vectorT), INTENT(IN) :: deps
+  TYPE(vectorT), INTENT(INOUT) :: dsig
+  TYPE(matrixT) :: Dmat   !> constitutive matrix
+
+  !> get constitutive matrix
+  CALL dmatrix(emod,nu, Dmat)
+
+  CALL vec_clean(dsig)    !> ensure dsig is deallocated (avoid memory leak)
+  dsig = Dmat*deps        !> compute stress increment
+
+  !> deallocate constitutive matrix
+  CALL dm_clean(Dmat)
+
+END SUBROUTINE linear_elastic_2d_plane_strain_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+! none
+
+END MODULE material_model

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/pde_solver_constants.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module defining PDE Solver Constants
+! ------------------------------------------------------------------------
+MODULE pde_solver_constants
+USE system_constants    !> Global system constants
+
+IMPLICIT NONE
+
+INTEGER, PARAMETER :: NGAUSS_ELEM = 1     !> number of Gaussian integration points per body element
+INTEGER, PARAMETER :: NGAUSS_BOUND = 1    !> number of Gaussian integration points per traction element
+
+!> Gaussian integration data for body elements
+DOUBLE PRECISION, PARAMETER, DIMENSION(NGAUSS_ELEM,2) :: GAUSS_PT_ELEM = RESHAPE( (/ ONE_THIRD, ONE_THIRD /), SHAPE(GAUSS_PT_ELEM) )
+DOUBLE PRECISION, PARAMETER, DIMENSION(NGAUSS_ELEM) :: GAUSS_WT_ELEM = (/ 1.d0 /)
+
+!> Gaussian integration data for traction elements
+DOUBLE PRECISION, PARAMETER, DIMENSION(NGAUSS_BOUND) :: GAUSS_PT_BOUND = (/ 0.5d0 /)
+DOUBLE PRECISION, PARAMETER, DIMENSION(NGAUSS_BOUND) :: GAUSS_WT_BOUND = (/ 1.d0 /)
+
+!> parameters for Newmark time-stepping
+DOUBLE PRECISION, PARAMETER :: GAMA = 0.5d0
+DOUBLE PRECISION, PARAMETER :: BETA = 0.25d0
+
+END MODULE pde_solver_constants

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/pde_solver_control.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for PDE Solver Control
+! ------------------------------------------------------------------------
+MODULE pde_solver_control
+USE system_constants              !> global system constants
+USE pde_solver_constants          !> constants for PDE solver
+USE log_message_control           !> print log/error messages
+USE log_messages                  !> log/error message and sender codes
+USE band_sym_matrix_def           !> banded symmetric matrix ADT
+USE dense_matrix_def              !> dense matrix ADT
+USE vector_def                    !> vector ADT
+USE field_data                    !> field data module
+USE boundary_data                 !> boundary data module
+USE material_data                 !> material property data module
+USE body_element_integration      !> body element integration (mass and stiff)
+USE traction_element_integration  !> traction element load vectors
+
+IMPLICIT NONE
+PRIVATE
+
+
+
+! ************************************************************************
+! EXPORTS
+! ************************************************************************
+
+
+!> Exported state variables
+PUBLIC :: hbw, nnod, nel, nelb, ndof,                 &
+          mass, modMass, damp, stiff,                 &
+          initStress, initStrain, body, trac, load,   &
+          prevDisp, incDisp, newDisp,                 &
+          prevVel, incVel, newVel,                    &
+          prevAcc, newAcc,                            &
+          prevStress, incStress, newStress,           &
+          prevStrain, incStrain, newStrain
+
+
+!> Exported interfaces
+PUBLIC :: pde_init, pde_clean, &
+          pde_buildMassMatrix, &
+          pde_buildStiffMatrix!, &
+!          pde_buildDampMatrix, &
+!          pde_buildModMassMatrix, &
+!          pde_buildLoadVector, &
+!          pde_initAcc, &
+!          pde_incAcc, pde_incVel, pde_incDisp, &
+!          pde_incStress, pde_incStrain, &
+!          pde_updateAcc, pde_updateVel, pde_updateDisp, &
+!          pde_updateStress, pde_updateStrain
+
+
+
+! ************************************************************************
+! LOCAL CONSTANTS
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! DATA TYPES
+! ************************************************************************
+
+! none
+
+
+
+! ************************************************************************
+! STATE VARIABLES
+! ************************************************************************
+
+INTEGER, SAVE               :: hbw            !> half bandwidth of system matrices
+INTEGER, SAVE               :: nnod           !> total number of nodes
+INTEGER, SAVE               :: nel            !> total number of body elements
+INTEGER, SAVE               :: nelb           !> total number of traction elements
+INTEGER, SAVE               :: ndof           !> total number of system degrees of freedom
+TYPE(bandSymMatrixT), SAVE  :: mass           !> mass matrix
+TYPE(bandSymMatrixT), SAVE  :: modMass        !> modified mass matrix (for time-stepping)
+TYPE(bandSymMatrixT), SAVE  :: stiff          !> stiffness matrix
+TYPE(bandSymMatrixT), SAVE  :: damp           !> damping matrix
+TYPE(vectorT), SAVE         :: initStress     !> load vector due to initial stress
+TYPE(vectorT), SAVE         :: initStrain     !> load vector due to initial strain
+TYPE(vectorT), SAVE         :: body           !> load vector due to body forces
+TYPE(vectorT), SAVE         :: trac           !> load vector due to surface tractions
+TYPE(vectorT), SAVE         :: load           !> total load vector
+TYPE(vectorT), SAVE         :: prevDisp       !> displacement vector on previous time step
+TYPE(vectorT), SAVE         :: incDisp        !> change in displacement on current time step
+TYPE(vectorT), SAVE         :: newDisp        !> displacement vector on current time step
+TYPE(vectorT), SAVE         :: prevVel        !> velocity vector on previous time step
+TYPE(vectorT), SAVE         :: incVel         !> change in velocity on current time step
+TYPE(vectorT), SAVE         :: newVel         !> velocity vector on current time step
+TYPE(vectorT), SAVE         :: prevAcc        !> acceleration vector on previous time step
+TYPE(vectorT), SAVE         :: newAcc         !> acceleration vector on current time step
+TYPE(vectorT), SAVE         :: prevStress     !> element stresses on previous time step
+TYPE(vectorT), SAVE         :: incStress      !> change in element stresses on current time step
+TYPE(vectorT), SAVE         :: newStress      !> element stresses on current time step
+TYPE(vectorT), SAVE         :: prevStrain     !> element strains on previous time step
+TYPE(vectorT), SAVE         :: incStrain      !> change in element strains on current time step
+TYPE(vectorT), SAVE         :: newStrain      !> element strains on current time step
+
+
+
+! ************************************************************************
+! INTERFACES
+! ************************************************************************
+
+!> \brief Interface to PDE solver initializer
+INTERFACE pde_init
+  MODULE PROCEDURE pde_init_
+END INTERFACE pde_init
+
+!> \brief Interface to PDE solver destructor
+INTERFACE pde_clean
+  MODULE PROCEDURE pde_clean_
+END INTERFACE pde_clean
+
+!> \brief Interface to mass matrix constructor
+INTERFACE pde_buildMassMatrix
+  MODULE PROCEDURE pde_build_mass_matrix_
+END INTERFACE pde_buildMassMatrix
+
+!> \brief Interface to stiffness matrix constructor
+INTERFACE pde_buildStiffMatrix
+  MODULE PROCEDURE pde_build_stiff_matrix_
+END INTERFACE pde_buildStiffMatrix
+
+!> \brief Interface to damping matrix constructor
+!INTERFACE pde_buildDampMatrix
+!  MODULE PROCEDURE pde_build_damp_matrix_
+!END INTERFACE pde_buildDampMatrix
+
+!> \brief Interface to modified mass matrix constructor
+!INTERFACE pde_buildModMassMatrix
+!  MODULE PROCEDURE pde_build_mod_mass_matrix_
+!END INTERFACE pde_buildModMassMatrix
+
+!> \brief Interface to load vector constructor
+!INTERFACE pde_buildLoadVector
+!  MODULE PROCEDURE pde_build_load_vector_
+!END INTERFACE pde_buildLoadVector
+
+!> \brief Interface to acceleration vector initializer
+!INTERFACE pde_initAcc
+!  MODULE PROCEDURE pde_init_acc_
+!END INTERFACE pde_initAcc
+
+!> \brief Interface to acceleration vector incrementer
+!INTERFACE pde_incAcc
+!  MODULE PROCEDURE pde_inc_acc_
+!END INTERFACE pde_incAcc
+
+!> \brief Interface to velocity vector incrementer
+!INTERFACE pde_incVel
+!  MODULE PROCEDURE pde_inc_vel_
+!END INTERFACE pde_incVel
+
+!> \brief Interface to displacement vector incrementer
+!INTERFACE pde_incDisp
+!  MODULE PROCEDURE pde_inc_disp_
+!END INTERFACE pde_incDisp
+
+!> \brief Interface to stress incrementer
+!INTERFACE pde_incStress
+!  MODULE PROCEDURE pde_inc_stress_
+!END INTERFACE pde_incStress
+
+!> \brief Interface to strain incrementer
+!INTERFACE pde_incStrain
+!  MODULE PROCEDURE pde_inc_strain_
+!END INTERFACE pde_incStrain
+
+!> \brief Interface to acceleration vector updater
+!INTERFACE pde_updateAcc
+!  MODULE PROCEDURE pde_update_acc_
+!END INTERFACE pde_updateAcc
+
+!> \brief Interface to velocity vector updater
+!INTERFACE pde_updateVel
+!  MODULE PROCEDURE pde_update_vel_
+!END INTERFACE pde_updateVel
+
+!> \brief Interface to displacement vector updater
+!INTERFACE pde_updateDisp
+!  MODULE PROCEDURE pde_update_disp_
+!END INTERFACE pde_updateDisp
+
+!> \brief Interface to stress updater
+!INTERFACE pde_updateStress
+!  MODULE PROCEDURE pde_update_stress_
+!END INTERFACE pde_updateStress
+
+!> \brief Interface to strain updater
+!INTERFACE pde_updateStrain
+!  MODULE PROCEDURE pde_update_strain_
+!END INTERFACE pde_updateStrain
+
+
+
+CONTAINS
+
+
+
+! ************************************************************************
+! ACCESS PROGRAMS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief Initializer for PDE solver
+!!
+!!  \param  exc For returning exception codes. (Optional).
+!!
+!!  This routine sets up the solution space for the problem by allocating
+!!  memory for the following state variables:
+!!
+!!    mass, stiff, initStress, initStrain, body, trac, load,
+!!    prevDisp, incDisp, newDisp,
+!!    prevVel, incVel, newVel,
+!!    prevAcc, newAcc,
+!!    prevStress, incStress, newStress,
+!!    prevStrain, incStrain, newStrain
+!!
+!!  Note that modMass and damp are not allocated here since
+!!  they will be computed by combining other matrices later.
+!!
+!!  This routine also initializes the values stored in:
+!!
+!!    prevDisp, prevVel, prevStress, prevStrain
+!!
+!!  These represent the initial conditions of the system.
+! ------------------------------------------------------------------------
+SUBROUTINE pde_init_ (exc)
+  INTEGER, INTENT(OUT), OPTIONAL :: exc
+  INTEGER :: exc_tmp    !> for storing exceptions
+  INTEGER :: nstress    !> total number of stress/strain entries
+  INTEGER :: i,j        !> loop variables
+  INTEGER :: dof        !> degree of freedom number
+
+  !> determine system size properties
+  nnod  = fld_numNode()
+  nel   = fld_numElem()
+  nelb  = bnd_numBoundElem()
+  ndof  = fld_numDof()
+  hbw   = compute_hbw()
+
+  ! ----------------------------
+  !> allocate system matrices
+  ! ----------------------------
+
+  !> mass matrix
+  CALL bsm_init(mass, hbw,ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> stiffness matrix
+  CALL bsm_init(stiff, hbw,ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+
+  ! ----------------------------
+  !> allocate system vectors
+  ! ----------------------------
+
+  !> initial stress load vector
+  CALL vec_init(initStress, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> initial strain load vector
+  CALL vec_init(initStrain, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> body force load vector
+  CALL vec_init(body, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> traction load vector
+  CALL vec_init(trac, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> total load vector
+  CALL vec_init(load, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> previous displacement
+  CALL vec_init(prevDisp, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> incremental displacement
+  CALL vec_init(incDisp, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> current displacement
+  CALL vec_init(newDisp, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> previous velocity
+  CALL vec_init(prevVel, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> incremental velocity
+  CALL vec_init(incVel, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> current velocity
+  CALL vec_init(newVel, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> previous acceleration
+  CALL vec_init(prevAcc, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> current acceleration
+  CALL vec_init(newAcc, ndof, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+
+  ! ----------------------------------
+  !> allocate stress/strain vectors
+  ! ----------------------------------
+
+  nstress = nel*NTNS
+
+  !> previous stress
+  CALL vec_init(prevStress, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> incremental stress
+  CALL vec_init(incStress, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> current stress
+  CALL vec_init(newStress, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> previous strain
+  CALL vec_init(prevStrain, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> incremental strain
+  CALL vec_init(incStrain, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+  !> current strain
+  CALL vec_init(newStrain, nstress, exc_tmp)
+  IF (exc_tmp.NE.OK) THEN
+    IF (PRESENT(exc)) exc=exc_tmp
+    RETURN
+  END IF
+
+
+  ! ----------------------------
+  !> set up initial conditions
+  ! ----------------------------
+
+  !> displacement and velocity
+  DO i = 1,nnod
+    DO j = 1,NDIM
+
+      dof = fld_getDof(i,j, exc_tmp)
+      IF (exc_tmp.NE.OK) THEN
+        IF (PRESENT(exc)) exc=exc_tmp
+        RETURN
+      END IF
+
+      IF (dof.NE.0) THEN
+
+        CALL vec_set(prevDisp, dof, fld_getDisp(i,j), exc_tmp)
+        IF (exc_tmp.NE.OK) THEN
+          IF (PRESENT(exc)) exc=exc_tmp
+          RETURN
+        END IF
+
+        CALL vec_set(prevVel, dof, fld_getVel(i,j), exc_tmp)
+        IF (exc_tmp.NE.OK) THEN
+          IF (PRESENT(exc)) exc=exc_tmp
+          RETURN
+        END IF
+
+      END IF
+
+    END DO
+  END DO
+
+  !> stress and strain
+  DO i = 1,nel
+
+    dof = (i-1)*NTNS
+
+    DO j = 1,NTNS
+
+      CALL vec_set(prevStress, dof+j, fld_getStressElem(i,j), exc_tmp)
+      IF (exc_tmp.NE.OK) THEN
+        IF (PRESENT(exc)) exc=exc_tmp
+        RETURN
+      END IF
+
+      CALL vec_set(prevStrain, dof+j, fld_getStrainElem(i,j), exc_tmp)
+      IF (exc_tmp.NE.OK) THEN
+        IF (PRESENT(exc)) exc=exc_tmp
+        RETURN
+      END IF
+
+    END DO
+
+  END DO
+
+
+  !> if no exceptions were raised, everything is A-OK
+  IF (PRESENT(exc)) exc=OK
+
+END SUBROUTINE pde_init_
+
+
+! ------------------------------------------------------------------------
+!>  \brief Destructor for PDE solver
+!!
+!!  This routine ensures that all state variables that have been
+!!  dynamically allocated are deallocated. It also resets system size
+!!  parameters.
+! ------------------------------------------------------------------------
+SUBROUTINE pde_clean_ ()
+
+  !> size parameters
+  nnod  = 0
+  nel   = 0
+  nelb  = 0
+  ndof  = 0
+  hbw   = 0
+
+  !> matrices
+  CALL bsm_clean(mass)
+  CALL bsm_clean(modMass)
+  CALL bsm_clean(stiff)
+  CALL bsm_clean(damp)
+
+  !> vectors
+  CALL vec_clean(initStress)
+  CALL vec_clean(initStrain)
+  CALL vec_clean(body)
+  CALL vec_clean(trac)
+  CALL vec_clean(load)
+  CALL vec_clean(prevDisp)
+  CALL vec_clean(incDisp)
+  CALL vec_clean(newDisp)
+  CALL vec_clean(prevVel)
+  CALL vec_clean(incVel)
+  CALL vec_clean(newVel)
+  CALL vec_clean(prevAcc)
+  CALL vec_clean(newAcc)
+  CALL vec_clean(prevStress)
+  CALL vec_clean(incStress)
+  CALL vec_clean(newStress)
+  CALL vec_clean(prevStrain)
+  CALL vec_clean(incStrain)
+  CALL vec_clean(newStrain)
+
+END SUBROUTINE pde_clean_
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Build the global mass matrix
+!!
+!!  \param  exc   Error code
+!!
+!!  This routine constructs the global (consistent) mass matrix by
+!!  summing element mass matrices taking connectivity into account. This
+!!  routine assumes that pde_init() has already been called.
+! ------------------------------------------------------------------------
+SUBROUTINE pde_build_mass_matrix_ (exc)
+  INTEGER, INTENT(OUT), OPTIONAL :: exc
+  TYPE(matrixT) :: emass      !> element mass matrix
+  INTEGER :: ind(NNODEL*NDIM) !> mapping indices
+  INTEGER :: i    !> loop variable
+  INTEGER :: exc_tmp    !> for getting exception code
+
+  !> loop through elements, adding element mass matrices
+  DO i = 1,nel
+
+    !> get element mass matrix and add it to the global mass matrix
+    CALL bint_emass(i,emass)
+    ind = get_index(i)
+    CALL bsm_mappedAdd(mass, emass,ind, exc_tmp)
+
+    !> check that mapped add was successful
+    IF (exc_tmp.NE.OK) THEN
+      IF (PRESENT(exc)) exc=exc_tmp
+      RETURN
+    END IF
+
+  END DO
+
+  !> if loop completed, everything is OK
+  IF (PRESENT(exc)) exc=OK
+
+END SUBROUTINE pde_build_mass_matrix_
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Build the global stiffness matrix
+!!
+!!  \param  exc   Error code
+!!
+!!  This routine constructs the global stiffness matrix by summing element
+!!  stiffness matrices taking connectivity into account. This
+!!  routine assumes that pde_init() has already been called.
+! ------------------------------------------------------------------------
+SUBROUTINE pde_build_stiff_matrix_ (exc)
+  INTEGER, INTENT(OUT), OPTIONAL :: exc
+  TYPE(matrixT) :: estiff      !> element stiffness matrix
+  INTEGER :: ind(NNODEL*NDIM) !> mapping indices
+  INTEGER :: i    !> loop variable
+  INTEGER :: exc_tmp    !> for getting exception code
+
+  !> loop through elements, adding element stiffness matrices
+  DO i = 1,nel
+
+    !> get element stiffness matrix and add it to the global stiffness matrix
+    CALL bint_estiff(i,estiff)
+    ind = get_index(i)
+    CALL bsm_mappedAdd(stiff, estiff,ind, exc_tmp)
+
+    !> check that mapped add was successful
+    IF (exc_tmp.NE.OK) THEN
+      IF (PRESENT(exc)) exc=exc_tmp
+      RETURN
+    END IF
+
+  END DO
+
+  !> if loop completed, everything is OK
+  IF (PRESENT(exc)) exc=OK
+
+END SUBROUTINE pde_build_stiff_matrix_
+
+
+
+! ************************************************************************
+! LOCAL FUNCTIONS
+! ************************************************************************
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Compute the half bandwidth of the system matrices
+!!
+!!  \return hbw   Half bandwidth
+! ------------------------------------------------------------------------
+FUNCTION compute_hbw () RESULT(hbw)
+  INTEGER :: hbw
+  INTEGER :: i,j,k,l,m  !> loop variables
+  INTEGER :: lmin, mmin !> loop limits
+  INTEGER :: nel        !> number of elements
+  INTEGER :: dof1,dof2  !> for calculating diff between dofs
+
+  !> initialize half bandwidth
+  hbw = 0
+
+  !> get number of elements
+  nel = fld_numElem()
+
+  !> loop through elements
+  !! (the goal is to find the maximum absolute difference between
+  !!  non-zero degree of freedom numbers within a single element)
+  DO i = 1,nel
+    DO j = 1,NNODEL   !> loop through nodes in current element
+      DO k = 1,NDIM   !> loop through coordinate (dof) directions
+
+        !> get current dof
+        dof1 = fld_getDof(fld_getConnect(i,j),k)
+
+        !> if current dof=0, no need to compare
+        IF (dof1.EQ.0)  CYCLE
+
+        !> get minimum loop indices for next dof
+        IF (k.LT.NDIM) THEN
+          lmin = j    !> start from current node
+          mmin = k+1  !> start from next dof of current node
+        ELSE
+          lmin = j+1  !> start from next node
+          mmin = 1    !> start from first dof of next node
+        END IF
+
+        DO l = lmin,NNODEL
+          DO m = mmin,NDIM
+
+            !> get next dof
+            dof2 = fld_getDof(fld_getConnect(i,l),m)
+
+            !> if next dof is non-zero, update hbw
+            IF (dof2.NE.0) THEN
+              hbw = MAX(hbw, ABS(dof2-dof1))
+            END IF
+
+          END DO  ! m = mmin,NDIM
+        END DO  ! l = lmin,NNODEL
+
+      END DO  ! k = 1,NDIM
+    END DO  ! j = 1,NNODEL
+  END DO  ! i = 1,nel
+
+  !> add one to include diagonal
+  hbw = hbw+1
+
+END FUNCTION compute_hbw
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Determine the mapping indices for a body element
+!!
+!!  \param  i     Element number
+!!
+!!  \return ind   Mapping indices
+!!
+!!  Note that this function assumes that the element index is valid. It
+!!  does not catch POSIT exceptions
+! ------------------------------------------------------------------------
+FUNCTION get_index (i) RESULT(ind)
+  INTEGER, INTENT(IN) :: i
+  INTEGER :: ind(NNODEL*NDIM)
+  INTEGER :: j,k    !> loop variables
+
+  DO j = 1,NNODEL
+    DO k = 1,NDIM
+      ind( (j-1)*NDIM+k ) = fld_getDof(fld_getConnect(i,j),k)
+    END DO
+  END DO
+
+END FUNCTION
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Determine the mapping indices for a traction element
+!!
+!!  \param  i     Element number
+!!
+!!  \return ind   Mapping indices
+!!
+!!  Note that this function assumes that the element index is valid. It
+!!  does not catch POSIT exceptions
+! ------------------------------------------------------------------------
+FUNCTION get_index_trac (i) RESULT(ind)
+  INTEGER, INTENT(IN) :: i
+  INTEGER :: ind(NNODELB*NDIM)
+  INTEGER :: j,k    !> loop variables
+
+  DO j = 1,NNODELB
+    DO k = 1,NDIM
+      ind( (j-1)*NDIM+k ) = fld_getDof(bnd_getConnect(i,j),k)
+    END DO
+  END DO
+
+END FUNCTION
+
+END MODULE pde_solver_control

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/pde_solver_control_test.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module for testing PDE Solver Control module
+! ------------------------------------------------------------------------
+MODULE pde_solver_control_test
+USE fruit                   !> FRUIT unit testing framework for Fortran
+USE system_constants        !> global system constants
+USE pde_solver_constants    !> constants for PDE solver
+USE log_message_control     !> print log/error messages
+USE log_messages            !> log/error message and sender codes
+USE band_sym_matrix_def     !> banded symmetric matrix ADT
+USE dense_matrix_def        !> dense matrix ADT
+USE vector_def              !> vector ADT
+USE field_data              !> field data module
+USE boundary_data           !> boundary data module
+USE material_data           !> material property data module
+USE pde_solver_control      !> PDE solver contorl module
+
+IMPLICIT NONE
+
+
+CONTAINS
+
+
+! ------------------------------------------------------------------------
+!>  \brief  Initialization test for PDE Solver Control
+!!
+!!  \param unit_name_size     Name of unit test (for FRUIT)
+!!  \param unit_name_area     Name of unit test (for FRUIT)
+!!  \param unit_name_mass     Name of unit test (for FRUIT)
+!!  \param unit_name_stiff    Name of unit test (for FRUIT)
+!!  \param unit_name_clean    Name of unit test (for FRUIT)
+!!  \param testName           Filename for log file (required for exceptions)
+!!  \param nnodExp            Expected number of nodes (15)
+!!  \param nelExp             Expected number of elements (16)
+!!  \param ndofExp            Expected number of degrees of freedom (24)
+!!  \param hbwExp             Expected half bandwidth (10)
+!!  \param ntnsExp            Expected number of tensor variables (48)
+!!  \param nmtl               Number of material types (1)
+!!  \param rho                Density of material 1 (2200 kg/m^3)
+!!  \param emod               Elastic modulus of material 1 (26000 MPa)
+!!  \param nu                 Poisson's ratio of material 1 (0.25)
+!!  \param expArea            Expected total area (5.0 x 8.0 = 40.0)
+!!  \param totalArea          Actual total area (from summing element areas)
+!!  \param expMass            For checking mass matrix (see test_straight.xlsx)
+!!  \param expStiff           For checking stiffness matrix (see test_straight.xlsx)
+!!  \param toler              Tolerance for floating point error
+!!  \param expMsg             Expected exception message
+!!  \param actMsg             Actual exception message
+!!
+!!  This routine tests tests correct allocation and initialization of
+!!  state variables for a simple system. The system is a straight column
+!!  with a fixed base. It has a width of 5 m and a height of 8 m. The
+!!  mesh of the system has the following properties:
+!!
+!!    nnod  = 15
+!!    nel   = 16
+!!    ndof  = 24
+!!    hbw   = 10
+!!
+!!  These parameters are verified after the initialization. The test also
+!!  verifies the size of all system matrices and vectors. The test also
+!!  deallocates the system and verifies that everything is set back to
+!!  zero post-deallocation.
+! ------------------------------------------------------------------------
+SUBROUTINE test_pde_solver_initialization
+  CHARACTER (LEN=*), PARAMETER :: unit_name_size = 'test_pde_solver_initialization_size'
+  CHARACTER (LEN=*), PARAMETER :: unit_name_area = 'test_pde_solver_initialization_area'
+  CHARACTER (LEN=*), PARAMETER :: unit_name_mass = 'test_pde_solver_initialization_mass'
+  CHARACTER (LEN=*), PARAMETER :: unit_name_stiff = 'test_pde_solver_initialization_stiff'
+  CHARACTER (LEN=*), PARAMETER :: unit_name_clean = 'test_pde_solver_initialization_clean'
+  CHARACTER (LEN=*), PARAMETER :: testName = 'test_straight_coarse'
+  INTEGER, PARAMETER :: nnodExp = 15
+  INTEGER, PARAMETER :: nelExp = 16
+  INTEGER, PARAMETER :: ndofExp = 24
+  INTEGER, PARAMETER :: hbwExp = 10
+  INTEGER, PARAMETER :: ntnsExp = nelExp*NTNS
+  INTEGER, PARAMETER :: nmtl = 1
+  DOUBLE PRECISION, PARAMETER :: rho = 2200.d0
+  DOUBLE PRECISION, PARAMETER :: emod = 26.d9
+  DOUBLE PRECISION, PARAMETER :: nu = 0.25d0
+  DOUBLE PRECISION, PARAMETER :: expArea = 40.d0
+  DOUBLE PRECISION :: totalArea
+  DOUBLE PRECISION, PARAMETER :: toler = 1.d-12
+  DOUBLE PRECISION :: expMass(hbwExp,ndofExp)
+  DOUBLE PRECISION :: expStiff(hbwExp,ndofExp)
+  INTEGER, PARAMETER :: expMsg=OK
+  INTEGER :: actMsg
+  INTEGER :: i,j    !> loop variables
+  DOUBLE PRECISION :: ycoord    !> for setting up coordinates
+
+  !> initialize log message file
+  CALL log_setFileName(testName)
+  CALL log_initLogFile()
+
+  !> initialize material data
+  CALL mtl_init(nmtl)
+  CALL mtl_setDens(1,rho)
+  CALL mtl_setEmod(1,emod)
+  CALL mtl_setPois(1,nu)
+
+  !> initialize node data
+  CALL fld_initNode(nnodExp)
+
+  !> x-coords of left column of nodes = 0
+
+  !> x-coords of middle column of nodes = 2.5
+  DO i = 2,nnodExp,3
+    CALL fld_setCoord(i,1,2.5d0)
+  END DO
+
+  !> x-coords of right column of nodes = 5
+  DO i = 3,nnodExp,3
+    CALL fld_setCoord(i,1,5.d0)
+  END DO
+
+  !> y-coords, steps of 2
+  DO i = 1,nnodExp,3
+    ycoord = ONE_THIRD*DBLE(i-1)*2.d0
+    CALL fld_setCoord(i,2,ycoord)
+    CALL fld_setCoord(i+1,2,ycoord)
+    CALL fld_setCoord(i+2,2,ycoord)
+  END DO
+
+  !> fix base
+  CALL fld_setFix(1,1,.TRUE.)
+  CALL fld_setFix(1,2,.TRUE.)
+  CALL fld_setFix(2,1,.TRUE.)
+  CALL fld_setFix(2,2,.TRUE.)
+  CALL fld_setFix(3,1,.TRUE.)
+  CALL fld_setFix(3,2,.TRUE.)
+
+  !> initialize system degrees of freedom
+  CALL fld_initDof()
+
+
+  !> initialize element data
+  CALL fld_initElem(nelExp)
+
+  !> element 1
+  CALL fld_setConnect(1,1,1)
+  CALL fld_setConnect(1,2,2)
+  CALL fld_setConnect(1,3,4)
+
+  !> element 2
+  CALL fld_setConnect(2,1,2)
+  CALL fld_setConnect(2,2,5)
+  CALL fld_setConnect(2,3,4)
+
+  !> element 3
+  CALL fld_setConnect(3,1,2)
+  CALL fld_setConnect(3,2,6)
+  CALL fld_setConnect(3,3,5)
+
+  !> element 4
+  CALL fld_setConnect(4,1,2)
+  CALL fld_setConnect(4,2,3)
+  CALL fld_setConnect(4,3,6)
+
+  !> element 5
+  CALL fld_setConnect(5,1,4)
+  CALL fld_setConnect(5,2,5)
+  CALL fld_setConnect(5,3,7)
+
+  !> element 6
+  CALL fld_setConnect(6,1,5)
+  CALL fld_setConnect(6,2,8)
+  CALL fld_setConnect(6,3,7)
+
+  !> element 7
+  CALL fld_setConnect(7,1,5)
+  CALL fld_setConnect(7,2,9)
+  CALL fld_setConnect(7,3,8)
+
+  !> element 8
+  CALL fld_setConnect(8,1,5)
+  CALL fld_setConnect(8,2,6)
+  CALL fld_setConnect(8,3,9)
+
+  !> element 9
+  CALL fld_setConnect(9,1,7)
+  CALL fld_setConnect(9,2,8)
+  CALL fld_setConnect(9,3,10)
+
+  !> element 10
+  CALL fld_setConnect(10,1,8)
+  CALL fld_setConnect(10,2,11)
+  CALL fld_setConnect(10,3,10)
+
+  !> element 11
+  CALL fld_setConnect(11,1,8)
+  CALL fld_setConnect(11,2,12)
+  CALL fld_setConnect(11,3,11)
+
+  !> element 12
+  CALL fld_setConnect(12,1,8)
+  CALL fld_setConnect(12,2,9)
+  CALL fld_setConnect(12,3,12)
+
+  !> element 13
+  CALL fld_setConnect(13,1,10)
+  CALL fld_setConnect(13,2,11)
+  CALL fld_setConnect(13,3,13)
+
+  !> element 14
+  CALL fld_setConnect(14,1,11)
+  CALL fld_setConnect(14,2,14)
+  CALL fld_setConnect(14,3,13)
+
+  !> element 15
+  CALL fld_setConnect(15,1,11)
+  CALL fld_setConnect(15,2,15)
+  CALL fld_setConnect(15,3,14)
+
+  !> element 16
+  CALL fld_setConnect(16,1,11)
+  CALL fld_setConnect(16,2,12)
+  CALL fld_setConnect(16,3,15)
+
+
+  !> set element materials
+  DO i = 1,nelExp
+    CALL fld_setMaterial(i,1)
+  END DO
+
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name_size)
+
+  !> initialize pde solver
+  CALL pde_init(actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> check system parameters
+  CALL assertEquals(nnodExp,nnod)
+  CALL assertEquals(nelExp,nel)
+  CALL assertEquals(ndofExp,ndof)
+  CALL assertEquals(hbwExp,hbw)
+
+  !> check state variables (matrices and vectors)
+  CALL assertEquals(hbwExp,bsm_halfBW(mass))
+  CALL assertEquals(ndofExp,bsm_numRows(mass))
+  CALL assertEquals(hbwExp,bsm_halfBW(stiff))
+  CALL assertEquals(ndofExp,bsm_numRows(stiff))
+  CALL assertEquals(ndofExp,vec_length(initStress))
+  CALL assertEquals(ndofExp,vec_length(initStrain))
+  CALL assertEquals(ndofExp,vec_length(body))
+  CALL assertEquals(ndofExp,vec_length(trac))
+  CALL assertEquals(ndofExp,vec_length(load))
+  CALL assertEquals(ndofExp,vec_length(prevDisp))
+  CALL assertEquals(ndofExp,vec_length(incDisp))
+  CALL assertEquals(ndofExp,vec_length(newDisp))
+  CALL assertEquals(ndofExp,vec_length(prevVel))
+  CALL assertEquals(ndofExp,vec_length(incVel))
+  CALL assertEquals(ndofExp,vec_length(newVel))
+  CALL assertEquals(ndofExp,vec_length(prevAcc))
+  CALL assertEquals(ndofExp,vec_length(newAcc))
+  CALL assertEquals(ntnsExp,vec_length(prevStress))
+  CALL assertEquals(ntnsExp,vec_length(incStress))
+  CALL assertEquals(ntnsExp,vec_length(newStress))
+  CALL assertEquals(ntnsExp,vec_length(prevStrain))
+  CALL assertEquals(ntnsExp,vec_length(incStrain))
+  CALL assertEquals(ntnsExp,vec_length(newStrain))
+
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name_area)
+
+  !> check total area
+  totalArea = 0.d0
+  DO i = 1,nel
+    totalArea = totalArea + fld_volElem(i)
+  END DO
+  CALL assertEquals(expArea,totalArea,toler)
+
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name_mass)
+
+  !> set up expected mass matrix data
+  expMass = RESHAPE( (/ &
+                      0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      3666.66666666666667d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1833.33333333333333d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      1222.22222222222222d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, 0.d0, 0.d0, &
+                      611.111111111111111d0, 0.d0, &
+                      1222.22222222222222d0 &
+                      /), SHAPE(expMass) )
+
+  !> set actual mass matrix
+  CALL pde_buildMassMatrix(actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> check mass matrix
+  CALL assertEquals(expMass,mass%dat, hbwExp,ndofExp, toler)
+
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name_stiff)
+
+  !> set up expected stiffness matrix data
+  expStiff = RESHAPE( (/ &
+                      0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 37960000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, 0.d0, 0.d0, 10400000000.d0, &
+                      47320000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, -24960000000.d0, -10400000000.d0, &
+                      75920000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, -10400000000.d0, -8320000000.d0, 0.d0, &
+                      94640000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 0.d0, -24960000000.d0, 10400000000.d0, &
+                      37960000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      0.d0, 10400000000.d0, -8320000000.d0, &
+                      -10400000000.d0, 47320000000.d0, 0.d0, 0.d0, &
+                      0.d0, -6500000000.d0, -5200000000.d0, 0.d0, &
+                      10400000000.d0, 0.d0, 0.d0, 37960000000.d0, &
+                      0.d0, 0.d0, -5200000000.d0, -19500000000.d0, &
+                      10400000000.d0, 0.d0, 0.d0, 0.d0, 10400000000.d0, &
+                      47320000000.d0, 0.d0, 0.d0, 0.d0, -13000000000.d0, &
+                      0.d0, 0.d0, 0.d0, -24960000000.d0, &
+                      -10400000000.d0, 75920000000.d0, 0.d0, 0.d0, 0.d0, &
+                      -39000000000.d0, 0.d0, 0.d0, -10400000000.d0, &
+                      -8320000000.d0, 0.d0, 94640000000.d0, 0.d0, 0.d0, &
+                      -10400000000.d0, -6500000000.d0, 5200000000.d0, &
+                      0.d0, 0.d0, -24960000000.d0, 10400000000.d0, &
+                      37960000000.d0, -10400000000.d0, 0.d0, &
+                      5200000000.d0, -19500000000.d0, 0.d0, 0.d0, &
+                      10400000000.d0, -8320000000.d0, -10400000000.d0, &
+                      47320000000.d0, 0.d0, 0.d0, 0.d0, -6500000000.d0, &
+                      -5200000000.d0, 0.d0, 10400000000.d0, 0.d0, 0.d0, &
+                      37960000000.d0, 0.d0, 0.d0, -5200000000.d0, &
+                      -19500000000.d0, 10400000000.d0, 0.d0, 0.d0, 0.d0, &
+                      10400000000.d0, 47320000000.d0, 0.d0, 0.d0, 0.d0, &
+                      -13000000000.d0, 0.d0, 0.d0, 0.d0, &
+                      -24960000000.d0, -10400000000.d0, 75920000000.d0, &
+                      0.d0, 0.d0, 0.d0, -39000000000.d0, 0.d0, 0.d0, &
+                      -10400000000.d0, -8320000000.d0, 0.d0, &
+                      94640000000.d0, 0.d0, 0.d0, -10400000000.d0, &
+                      -6500000000.d0, 5200000000.d0, 0.d0, 0.d0, &
+                      -24960000000.d0, 10400000000.d0, 37960000000.d0, &
+                      -10400000000.d0, 0.d0, 5200000000.d0, &
+                      -19500000000.d0, 0.d0, 0.d0, 10400000000.d0, &
+                      -8320000000.d0, -10400000000.d0, 47320000000.d0, &
+                      0.d0, 0.d0, 0.d0, -6500000000.d0, -5200000000.d0, &
+                      0.d0, 10400000000.d0, 0.d0, 0.d0, 18980000000.d0, &
+                      0.d0, 0.d0, -5200000000.d0, -19500000000.d0, &
+                      10400000000.d0, 0.d0, 0.d0, 0.d0, 0.d0, &
+                      23660000000.d0, 0.d0, 0.d0, 0.d0, -13000000000.d0, &
+                      0.d0, 0.d0, 0.d0, -12480000000.d0, -5200000000.d0, &
+                      37960000000.d0, 0.d0, 0.d0, 0.d0, -39000000000.d0, &
+                      0.d0, 0.d0, -5200000000.d0, -4160000000.d0, 0.d0, &
+                      47320000000.d0, 0.d0, 0.d0, -10400000000.d0, &
+                      -6500000000.d0, 5200000000.d0, 0.d0, 0.d0, &
+                      -12480000000.d0, 5200000000.d0, 18980000000.d0, &
+                      -10400000000.d0, 0.d0, 5200000000.d0, &
+                      -19500000000.d0, 0.d0, 0.d0, 5200000000.d0, &
+                      -4160000000.d0, 0.d0, 23660000000.d0 &
+                      /), SHAPE(expStiff) )
+
+  !> set actual stiffness matrix
+  CALL pde_buildStiffMatrix(actMsg)
+  CALL assertEquals(expMsg,actMsg)
+
+  !> check stiffness matrix
+  CALL assertEquals(expStiff,stiff%dat, hbwExp,ndofExp, toler)
+
+
+  !> initialize unit test
+  CALL set_unit_name(unit_name_clean)
+
+  !> finalize pde solver
+  CALL pde_clean()
+
+  !> check system parameters
+  CALL assertEquals(0,nnod)
+  CALL assertEquals(0,nel)
+  CALL assertEquals(0,ndof)
+  CALL assertEquals(0,hbw)
+
+  !> check state variables (matrices and vectors)
+  CALL assertEquals(0,bsm_halfBW(mass))
+  CALL assertEquals(0,bsm_numRows(mass))
+  CALL assertEquals(0,bsm_halfBW(stiff))
+  CALL assertEquals(0,bsm_numRows(stiff))
+  CALL assertEquals(0,vec_length(initStress))
+  CALL assertEquals(0,vec_length(initStrain))
+  CALL assertEquals(0,vec_length(body))
+  CALL assertEquals(0,vec_length(trac))
+  CALL assertEquals(0,vec_length(load))
+  CALL assertEquals(0,vec_length(prevDisp))
+  CALL assertEquals(0,vec_length(incDisp))
+  CALL assertEquals(0,vec_length(newDisp))
+  CALL assertEquals(0,vec_length(prevVel))
+  CALL assertEquals(0,vec_length(incVel))
+  CALL assertEquals(0,vec_length(newVel))
+  CALL assertEquals(0,vec_length(prevAcc))
+  CALL assertEquals(0,vec_length(newAcc))
+  CALL assertEquals(0,vec_length(prevStress))
+  CALL assertEquals(0,vec_length(incStress))
+  CALL assertEquals(0,vec_length(newStress))
+  CALL assertEquals(0,vec_length(prevStrain))
+  CALL assertEquals(0,vec_length(incStrain))
+  CALL assertEquals(0,vec_length(newStrain))
+
+  !> finalize log file and deallocate objects
+  CALL log_closeLogFile()
+  CALL fld_cleanNode()
+  CALL fld_cleanElem()
+
+END SUBROUTINE test_pde_solver_initialization
+
+END MODULE pde_solver_control_test

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/system_constants.f90.svn-base

+! ------------------------------------------------------------------------
+!>  \brief Module defining System Constants
+! ------------------------------------------------------------------------
+MODULE system_constants
+
+IMPLICIT NONE
+
+INTEGER, PARAMETER :: MAXLEN = 200            !> maximum string length
+
+DOUBLE PRECISION, PARAMETER :: ONE_THIRD = 0.33333333333333333d0    !> avoid computing 1/3
+
+INTEGER, PARAMETER :: NDIM          = 2       !> number of coordinate dims
+INTEGER, PARAMETER :: NTNS          = 3       !> number of tensor components
+INTEGER, PARAMETER :: NNODEL        = 3       !> number of nodes per body element
+INTEGER, PARAMETER :: NNODELB       = 2       !> number of nodes per traction element
+INTEGER, PARAMETER :: MAX_NODES     = 2000    !> maximum number of nodes
+INTEGER, PARAMETER :: MAX_DOFS      = 3990    !> maximum number of degrees of freedom
+INTEGER, PARAMETER :: MAX_ELEMENTS  = 5000    !> maximum number of elements
+INTEGER, PARAMETER :: MAX_BOUNDELS  = 2000    !> maximum number of traction elements
+INTEGER, PARAMETER :: MAX_MATERIALS = 30      !> maximum number of materials
+INTEGER, PARAMETER :: MAX_TIMESTEPS = 10000   !> maximum number of time steps
+
+DOUBLE PRECISION, PARAMETER :: E_MIN      =  0.d0       !> minimum value of elastic modulus
+DOUBLE PRECISION, PARAMETER :: E_MAX      =  1.d11      !> maximum value of elastic modulus
+DOUBLE PRECISION, PARAMETER :: NU_MIN     =  0.d0       !> minimum value of Poisson's ratio
+DOUBLE PRECISION, PARAMETER :: NU_MAX     =  0.499d0    !> maximum value of Poisson's ratio
+DOUBLE PRECISION, PARAMETER :: RHO_MIN    =  0.d0       !> minimum value of density
+DOUBLE PRECISION, PARAMETER :: RHO_MAX    =  1.d11      !> maximum value of density
+DOUBLE PRECISION, PARAMETER :: COORD_MIN  = -1.d11      !> minimum value of coordinates
+DOUBLE PRECISION, PARAMETER :: COORD_MAX  =  1.d11      !> maximum value of coordinates
+DOUBLE PRECISION, PARAMETER :: DISP_MIN   = -1.d11      !> minimum value of displacement
+DOUBLE PRECISION, PARAMETER :: DISP_MAX   =  1.d11      !> maximum value of displacement
+DOUBLE PRECISION, PARAMETER :: VEL_MIN    = -1.d11      !> minimum value of velocity
+DOUBLE PRECISION, PARAMETER :: VEL_MAX    =  1.d11      !> maximum value of velocity
+DOUBLE PRECISION, PARAMETER :: ACC_MIN    = -1.d11      !> minimum value of acceleration
+DOUBLE PRECISION, PARAMETER :: ACC_MAX    =  1.d11      !> maximum value of acceleration
+DOUBLE PRECISION, PARAMETER :: SIG_MIN    = -1.d11      !> minimum value of stress
+DOUBLE PRECISION, PARAMETER :: SIG_MAX    =  1.d11      !> maximum value of stress
+DOUBLE PRECISION, PARAMETER :: STR_SMALL  =  1.d-2      !> threshold for small strains
+DOUBLE PRECISION, PARAMETER :: STR_MIN    = -STR_SMALL  !> minimum value of strain
+DOUBLE PRECISION, PARAMETER :: STR_MAX    =  STR_SMALL  !> maximum value of strain
+DOUBLE PRECISION, PARAMETER :: DTIME_MIN  =  1.d-11     !> minimum value of time step
+DOUBLE PRECISION, PARAMETER :: DTIME_MAX  =  1.d4       !> maximum value of time step
+
+END MODULE system_constants

Examples/SoilWaterStructureInteraction/dynsws/trunk/.svn/text-base/test_driver.f90.svn-base

+PROGRAM test_driver
+USE fruit
+USE log_messages_test
+USE log_message_control_test
+USE dense_matrix_test
+USE vector_test
+USE band_sym_matrix_test
+USE linear_solver_test
+USE material_data_test
+USE field_data_test
+USE boundary_data_test
+USE constitutive_test
+USE pde_solver_control_test
+
+
+  CALL init_fruit
+
+  ! log_messages message code tests
+  CALL test_ALLOC_message_code
+  CALL test_DIMEN_message_code
+  CALL test_EXCEED_message_code
+  CALL test_EXISTS_message_code
+  CALL test_FORMT_message_code
+  CALL test_POSDEF_message_code
+  CALL test_POSIT_message_code
+  CALL test_SZE_message_code
+  CALL test_TYP_message_code
+  CALL test_unexpected_message_code
+
+  ! log_messages sender code tests
+  CALL test_BFCRDR_sender_code
+  CALL test_BNDDAT_sender_code
+  CALL test_BNDRDR_sender_code
+  CALL test_BSYMAT_sender_code
+  CALL test_CNSMAT_sender_code
+  CALL test_DMNRDR_sender_code
+  CALL test_DNSMAT_sender_code
+  CALL test_FLDDAT_sender_code
+  CALL test_ICTRDR_sender_code
+  CALL test_ICVRDR_sender_code
+  CALL test_KBCRDR_sender_code
+  CALL test_LINSLV_sender_code
+  CALL test_MTLDAT_sender_code
+  CALL test_MTLRDR_sender_code
+  CALL test_NBCRDR_sender_code
+  CALL test_TNSWTR_sender_code
+  CALL test_VECTOR_sender_code
+  CALL test_VECWTR_sender_code
+  CALL test_unexpected_sender_code
+
+  ! log_message_control tests
+  CALL test_log_setFileName
+  CALL test_log_initLogFile
+  CALL test_log_closeLogFile
+  CALL test_log_printLogMsg
+
+  ! dense_matrix_def tests
+  CALL test_dm_allocation_MSG
+  CALL test_dm_allocation_SZE
+  CALL test_dm_allocation_DAT
+  CALL test_dm_deallocation
+  CALL test_dm_num_rows_not_allocated
+  CALL test_dm_num_rows_allocated
+  CALL test_dm_num_cols_not_allocated
+  CALL test_dm_num_cols_allocated
+  CALL test_dm_get_POSIT
+  CALL test_dm_get_VAL
+  CALL test_dm_set_POSIT
+  CALL test_dm_set_VAL
+  CALL test_dm_add_DIMEN
+  CALL test_dm_add_VAL
+  CALL test_dm_add_OP
+  CALL test_dm_scal_mul_VAL
+  CALL test_dm_scal_mul_ZERO
+  CALL test_dm_scal_mul_OP
+  CALL test_dm_vec_mul_DIMEN
+  CALL test_dm_vec_mul_VAL
+  CALL test_dm_vec_mul_ZERO
+  CALL test_dm_vec_mul_IDENT
+  CALL test_dm_vec_mul_OP
+  CALL test_dm_mat_mul_DIMEN
+  CALL test_dm_mat_mul_VAL
+  CALL test_dm_mat_mul_ZERO
+  CALL test_dm_mat_mul_IDENT
+  CALL test_dm_mat_mul_OP
+  CALL test_dm_transpose
+
+  ! vector_def tests
+  CALL test_vec_allocation_MSG
+  CALL test_vec_allocation_SZE
+  CALL test_vec_allocation_DAT
+  CALL test_vec_deallocation
+  CALL test_vec_length_not_allocated
+  CALL test_vec_length_allocated
+  CALL test_vec_get_POSIT
+  CALL test_vec_get_VAL
+  CALL test_vec_set_POSIT
+  CALL test_vec_set_VAL
+  CALL test_vec_add_DIMEN
+  CALL test_vec_add_VAL
+  CALL test_vec_add_OP
+  CALL test_vec_mapped_add_DIMEN
+  CALL test_vec_mapped_add_POSIT
+  CALL test_vec_mapped_add_VAL1
+  CALL test_vec_mapped_add_VAL2
+  CALL test_vec_mapped_add_VAL3
+  CALL test_vec_scal_mul_VAL
+  CALL test_vec_scal_mul_ZERO
+  CALL test_vec_scal_mul_OP
+  CALL test_vec_dot_prod_DIMEN
+  CALL test_vec_dot_prod_ZERO
+  CALL test_vec_dot_prod_VAL
+
+  ! band_sym_matrix_def tests
+  CALL test_bsm_allocation_MSG
+  CALL test_bsm_allocation_SZE
+  CALL test_bsm_allocation_DAT
+  CALL test_bsm_deallocation
+  CALL test_bsm_num_rows_not_allocated
+  CALL test_bsm_num_rows_allocated
+  CALL test_bsm_half_bw_not_allocated
+  CALL test_bsm_half_bw_allocated
+  CALL test_bsm_get_POSIT
+  CALL test_bsm_get_VAL
+  CALL test_bsm_set_POSIT
+  CALL test_bsm_set_VAL
+  CALL test_bsm_set_decomp_DIMEN
+  CALL test_bsm_set_decomp_VAL
+  CALL test_bsm_is_decomposed
+  CALL test_bsm_add_DIMEN
+  CALL test_bsm_add_VAL
+  CALL test_bsm_add_OP
+  CALL test_bsm_mapped_add_DIMEN1
+  CALL test_bsm_mapped_add_DIMEN2
+  CALL test_bsm_mapped_add_DIMEN3
+  CALL test_bsm_mapped_add_POSIT1
+  CALL test_bsm_mapped_add_POSIT2
+  CALL test_bsm_mapped_add_VAL1
+  CALL test_bsm_mapped_add_VAL2
+  CALL test_bsm_mapped_add_VAL3
+  CALL test_bsm_scal_mul_VAL
+  CALL test_bsm_scal_mul_ZERO
+  CALL test_bsm_scal_mul_OP
+  CALL test_bsm_vec_mul_DIMEN
+  CALL test_bsm_vec_mul_VAL
+  CALL test_bsm_vec_mul_ZERO
+  CALL test_bsm_vec_mul_IDENT
+  CALL test_bsm_vec_mul_OP
+
+  ! linear_solver tests
+  CALL test_linear_solver_DIMEN
+  CALL test_linear_solver_POSDEF
+  CALL test_linear_solver_VAL
+  CALL test_linear_solver_IDENT
+  CALL test_linear_solver_ZERO
+  CALL test_linear_solver_DECOMP
+
+  ! material property data tests
+  CALL test_mtl_allocation_MSG
+  CALL test_mtl_allocation_SZE
+  CALL test_mtl_num_mtl_not_allocated
+  CALL test_mtl_num_mtl_allocated
+  CALL test_mtl_get_emod_POSIT
+  CALL test_mtl_set_emod_POSIT
+  CALL test_mtl_set_emod_EXCEED
+  CALL test_mtl_get_set_emod_VAL
+  CALL test_mtl_get_pois_POSIT
+  CALL test_mtl_set_pois_POSIT
+  CALL test_mtl_set_pois_EXCEED
+  CALL test_mtl_get_set_pois_VAL
+  CALL test_mtl_get_dens_POSIT
+  CALL test_mtl_set_dens_POSIT
+  CALL test_mtl_set_dens_EXCEED
+  CALL test_mtl_get_set_dens_VAL
+
+  ! field data tests
+  CALL test_fld_init_time_EXCEED
+  CALL test_fld_time_step_not_initialized
+  CALL test_fld_time_step_initialized
+  CALL test_fld_num_time_step_not_initialized
+  CALL test_fld_num_time_step_initialized
+  CALL test_fld_node_allocation_MSG
+  CALL test_fld_node_allocation_SZE
+  CALL test_fld_num_node_not_allocated
+  CALL test_fld_num_node_allocated
+  CALL test_fld_get_coord_POSIT
+  CALL test_fld_set_coord_POSIT
+  CALL test_fld_set_coord_EXCEED
+  CALL test_fld_get_set_coord_VAL
+  CALL test_fld_get_fix_POSIT
+  CALL test_fld_set_fix_POSIT
+  CALL test_fld_get_set_fix_VAL
+  CALL test_fld_dof_initialization_SZE
+  CALL test_fld_num_dof_not_initialized
+  CALL test_fld_num_dof_initialized
+  CALL test_fld_get_dof_POSIT
+  CALL test_fld_get_dof_VAL
+  CALL test_fld_get_disp_POSIT
+  CALL test_fld_set_disp_POSIT
+  CALL test_fld_set_disp_EXCEED
+  CALL test_fld_get_set_disp_VAL
+  CALL test_fld_get_vel_POSIT
+  CALL test_fld_set_vel_POSIT
+  CALL test_fld_set_vel_EXCEED
+  CALL test_fld_get_set_vel_VAL
+  CALL test_fld_get_acc_POSIT
+  CALL test_fld_set_acc_POSIT
+  CALL test_fld_set_acc_EXCEED
+  CALL test_fld_get_set_acc_VAL
+  CALL test_fld_get_body_acc_POSIT
+  CALL test_fld_set_body_acc_POSIT
+  CALL test_fld_set_body_acc_EXCEED
+  CALL test_fld_get_set_body_acc_VAL
+  CALL test_fld_get_stress_node_POSIT
+  CALL test_fld_set_stress_node_POSIT
+  CALL test_fld_set_stress_node_EXCEED
+  CALL test_fld_get_set_stress_node_VAL
+  CALL test_fld_get_strain_node_POSIT
+  CALL test_fld_set_strain_node_POSIT
+  CALL test_fld_set_strain_node_EXCEED
+  CALL test_fld_get_set_strain_node_VAL
+  CALL test_fld_elem_allocation_MSG
+  CALL test_fld_elem_allocation_SZE
+  CALL test_fld_num_elem_not_allocated
+  CALL test_fld_num_elem_allocated
+  CALL test_fld_get_connect_POSIT
+  CALL test_fld_set_connect_POSIT
+  CALL test_fld_set_connect_EXCEED
+  CALL test_fld_get_set_connect_VAL
+  CALL test_fld_vol_elem_POSIT
+  CALL test_fld_vol_elem_ZERO
+  CALL test_fld_vol_elem_VAL
+  CALL test_fld_get_stress_elem_POSIT
+  CALL test_fld_set_stress_elem_POSIT
+  CALL test_fld_set_stress_elem_EXCEED
+  CALL test_fld_get_set_stress_elem_VAL
+  CALL test_fld_get_strain_elem_POSIT
+  CALL test_fld_set_strain_elem_POSIT
+  CALL test_fld_set_strain_elem_EXCEED
+  CALL test_fld_get_set_strain_elem_VAL
+
+  ! boundary data tests
+  CALL test_bnd_elem_allocation_MSG
+  CALL test_bnd_elem_allocation_SZE
+  CALL test_bnd_num_elem_not_allocated
+  CALL test_bnd_num_elem_allocated
+  CALL test_bnd_get_connect_POSIT
+  CALL test_bnd_set_connect_POSIT
+  CALL test_bnd_set_connect_EXCEED
+  CALL test_bnd_get_set_connect_VAL
+  CALL test_bnd_len_bound_elem_POSIT
+  CALL test_bnd_len_bound_elem_ZERO
+  CALL test_bnd_len_bound_elem_VAL
+  CALL test_bnd_get_trac_POSIT
+  CALL test_bnd_set_trac_POSIT
+  CALL test_bnd_set_trac_EXCEED
+  CALL test_bnd_get_set_trac_VAL
+
+  ! constitutive matrix tests
+  CALL test_constitutive_EXCEED
+  CALL test_constitutive_ZERO
+  CALL test_constitutive_VAL
+
+  ! pde solver tests
+  CALL test_pde_solver_initialization
+
+  CALL fruit_summary
+
+END PROGRAM test_driver