H.1.5 Add_to_Matrix_Equation_Monomial Procedure

The main documentation of the Add_to_Matrix_Equation_Monomial Procedure contains additional explanation of this code listing.

  subroutine Add_to_Matrix_Equation_Monomial (Monomial, ELLM, RHS_MV)

    ! Use association information.

    use Caesar_Numbers_Module, only: one, zero

    ! Input variables.

    type(Monomial_type), intent(inout) :: Monomial  ! Monomial to be added.

    ! Input/Output variables.

    ! ELL_Matrix to be incremented.
    type(ELL_Matrix_type), intent(inout) :: ELLM
    ! RHS mathematic vector.
    type(Mathematic_Vector_type), intent(inout) :: RHS_MV   

    ! Internal variables.

    ! Arrays to manipulate values of the matrix.
    type(real,2) :: Matrix_Values   
    type(integer,1) :: Matrix_Rows
    type(integer,2) :: Matrix_Columns
    type(integer) :: i                          ! Loop parameter.
    type(integer) :: shift                      ! Index shift.
    type(real,1) :: Volume_Cells   ! Volume of the cells.

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Verify requirements.

    VERIFY(Valid_State(Monomial),5) ! Monomial is valid.
    VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
    VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

    ! Create some working vectors.

    call Initialize (Matrix_Values, NRows_PE(ELLM), 1)
    call Initialize (Matrix_Rows, NRows_PE(ELLM))
    call Initialize (Matrix_Columns, NRows_PE(ELLM), 1)
    call Initialize (Volume_Cells, NRows_PE(ELLM))

    ! Increment ELLM values.

    call Get_Volume_Cells (Volume_Cells, Monomial%Mesh)
    if (Monomial%Exponent /= zero) then
      Matrix_Values(:,1) = Monomial%Coefficient * Monomial%Exponent * &
                           Monomial%Phi ** (Monomial%Exponent - one) * &
                           Volume_Cells
    else
      ! Treat exponent=0 (monomial=constant) case separately to avoid 
      ! division by zero if Phi=0 and to save a little time.
      Matrix_Values(:,1) = zero
    end if
    shift = First_PE(Monomial%Structure) - 1
    do i = 1, Length_PE(Monomial%Structure)
      Matrix_Rows(i) = i
      Matrix_Columns(i,1) = i + shift
    end do
    call Add_Values (ELLM, Matrix_Values, Matrix_Rows, Matrix_Columns, &
                     Global=.false.)

    ! Increment RHS_MV values.

    Matrix_Values(:,1) = Monomial%Coefficient * &
                         (Monomial%Exponent - one) * &
                         Monomial%Phi ** Monomial%Exponent * &
                         Volume_Cells
    call Add_Values (RHS_MV, Matrix_Values(:,1))

    ! Finalize working vectors.

    call Finalize (Matrix_Values)
    call Finalize (Matrix_Rows)
    call Finalize (Matrix_Columns)
    call Finalize (Volume_Cells)

    ! Verify guarantees.

    VERIFY(Valid_State(ELLM),5)     ! ELLM is valid.
    VERIFY(Valid_State(RHS_MV),5)   ! RHS_MV is valid.

    return
  end subroutine Add_to_Matrix_Equation_Monomial