Poiseuille_equations_routines.f90

Go to the documentation of this file.

MODULE poiseuille_equations_routines

  USE base_routines 
  USE basis_routines
  USE boundary_conditions_routines
  USE constants
  USE control_loop_routines
  USE distributed_matrix_vector
  USE domain_mappings
  USE equations_routines
  USE equations_mapping_routines
  USE equations_matrices_routines
  USE equations_set_constants
  USE field_routines
  USE input_output
  USE iso_varying_string
  USE kinds
  USE matrix_vector
  USE problem_constants
  USE strings
  USE solver_routines
  USE timer
  USE types
! temporary input for vector-source
  USE fluid_mechanics_io_routines

#include "macros.h"  

  IMPLICIT NONE

  PRIVATE

  !Module parameters

  !Module types

  !Module variables

  !Interfaces

  PUBLIC poiseuille_equation_equations_set_setup

  PUBLIC poiseuille_equationssetsolutionmethodset

  PUBLIC poiseuille_boundaryconditionsanalyticcalculate

  PUBLIC poiseuille_equationssetspecificationset

  PUBLIC poiseuille_finiteelementcalculate

  PUBLIC poiseuille_equation_problem_setup

  PUBLIC poiseuille_pre_solve,poiseuille_post_solve

  PUBLIC poiseuille_problemspecificationset
  
CONTAINS

  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_boundaryconditionsanalyticcalculate(EQUATIONS_SET,BOUNDARY_CONDITIONS,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(boundary_conditions_type), POINTER :: BOUNDARY_CONDITIONS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,deriv_idx,dim_idx,local_ny,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type
    INTEGER(INTG) :: COUNT_DOF
    REAL(DP) :: VALUE,X(3)
    REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:)
    TYPE(domain_type), POINTER :: DOMAIN
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE
    TYPE(varying_string) :: LOCAL_ERROR    
    
    enters("Poiseuille_BoundaryConditionsAnalyticCalculate",err,error,*999)
    
    IF(ASSOCIATED(equations_set)) THEN
      IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
        dependent_field=>equations_set%DEPENDENT%DEPENDENT_FIELD
        IF(ASSOCIATED(dependent_field)) THEN
          geometric_field=>equations_set%GEOMETRY%GEOMETRIC_FIELD
          IF(ASSOCIATED(geometric_field)) THEN            
            CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
            NULLIFY(geometric_variable)
            CALL field_variable_get(geometric_field,field_u_variable_type,geometric_variable,err,error,*999)
            CALL field_parameter_set_data_get(geometric_field,field_u_variable_type,field_values_set_type,geometric_parameters, &
              & err,error,*999)

            IF(ASSOCIATED(boundary_conditions)) THEN

              count_dof=0

              DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES !U and deludeln
                variable_type=dependent_field%VARIABLES(variable_idx)%VARIABLE_TYPE
                field_variable=>dependent_field%VARIABLE_TYPE_MAP(variable_type)%PTR
                IF(ASSOCIATED(field_variable)) THEN
                  CALL field_parameter_set_create(dependent_field,variable_type,field_analytic_values_set_type,err,error,*999)
                  DO component_idx=1,field_variable%NUMBER_OF_COMPONENTS !u,v,w
                    IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
                      domain=>field_variable%COMPONENTS(component_idx)%DOMAIN
                      IF(ASSOCIATED(domain)) THEN
                        IF(ASSOCIATED(domain%TOPOLOGY)) THEN
                          domain_nodes=>domain%TOPOLOGY%NODES
                          IF(ASSOCIATED(domain_nodes)) THEN

                            !Loop over the local nodes excluding the ghosts.
                            DO node_idx=1,domain_nodes%NUMBER_OF_NODES
                              !!TODO \todo We should interpolate the geometric field here and the node position.
                              DO dim_idx=1,number_of_dimensions
                                !Default to version 1 of each node derivative
                                local_ny=geometric_variable%COMPONENTS(dim_idx)%PARAM_TO_DOF_MAP%NODE_PARAM2DOF_MAP% &
                                  & nodes(node_idx)%DERIVATIVES(1)%VERSIONS(1)
                                x(dim_idx)=geometric_parameters(local_ny)
                              ENDDO !dim_idx
                              !Loop over the derivatives
                              DO deriv_idx=1,domain_nodes%NODES(node_idx)%NUMBER_OF_DERIVATIVES
                                SELECT CASE(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE)
                                CASE(equations_set_poiseuille_equation_two_dim_1)
                                  !u=ln(4/(x+y+1^2))
                                  SELECT CASE(variable_type)
                                  CASE(field_u_variable_type)
                                    SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
                                    CASE(no_global_deriv)
                                      VALUE=log(4.0_dp/((x(1)+x(2)+1.0_dp)**2))
                                    CASE(global_deriv_s1)
                                      CALL flagerror("Not implemented.",err,error,*999)
                                    CASE(global_deriv_s2)
                                      CALL flagerror("Not implemented.",err,error,*999)
                                    CASE(global_deriv_s1_s2)
                                      CALL flagerror("Not implemented.",err,error,*999)
                                    CASE DEFAULT
                                      local_error="The global derivative index of "//trim(number_to_vstring( &
                                      domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
                                        & err,error))//" is invalid."
                                      CALL flagerror(local_error,err,error,*999)
                                    END SELECT
                                  CASE(field_deludeln_variable_type)
                                    SELECT CASE(domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX)
                                    CASE(no_global_deriv)
                                      !This is believed to be incorrect, should be: VALUE=-2.0_DP/(X(1)+X(2)+1.0_DP)
                                      VALUE=-2.0_dp*(x(1)+x(2))/(x(1)+x(2)+1.0_dp)
                                    CASE(global_deriv_s1)
                                      CALL flagerror("Not implemented.",err,error,*999)
                                    CASE(global_deriv_s2)
                                      CALL flagerror("Not implemented.",err,error,*999)
                                    CASE(global_deriv_s1_s2)
                                     CALL flagerror("Not implemented.",err,error,*999)
                                    CASE DEFAULT
                                      local_error="The global derivative index of "//trim(number_to_vstring( &
                                      domain_nodes%NODES(node_idx)%DERIVATIVES(deriv_idx)%GLOBAL_DERIVATIVE_INDEX,"*", &
                                        & err,error))//" is invalid."
                                      CALL flagerror(local_error,err,error,*999)
                                    END SELECT
                                  END SELECT
                                CASE DEFAULT
                                  IF(variable_type==field_u_variable_type.AND.domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN !For Dirichlet
                                    CALL boundary_conditions_set_local_dof(boundary_conditions,dependent_field,variable_type, &
                                      & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
                                  ENDIF
                                END SELECT

                                IF(variable_type==field_deludeln_variable_type.and.node_idx/=1) THEN
                                  IF(domain_nodes%NODES(node_idx)%BOUNDARY_NODE) THEN
                                    !If we are a boundary node then set the analytic value on the boundary
                                    !CALL BOUNDARY_CONDITIONS_SET_LOCAL_DOF(BOUNDARY_CONDITIONS,variable_type,local_ny, &
                                    !  & BOUNDARY_CONDITION_FIXED,VALUE,ERR,ERROR,*999)
                                    !Do nothing at present
                                  ENDIF
                                ENDIF

                              ENDDO !deriv_idx
                            ENDDO !node_idx
                          ELSE
                            CALL flagerror("Domain topology nodes is not associated.",err,error,*999)
                          ENDIF
                        ELSE
                          CALL flagerror("Domain topology is not associated.",err,error,*999)
                        ENDIF
                      ELSE
                        CALL flagerror("Domain is not associated.",err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Only node based interpolation is implemented.",err,error,*999)
                    ENDIF
                  ENDDO !component_idx
                  CALL field_parameter_set_update_start(dependent_field,variable_type,field_analytic_values_set_type, &
                    & err,error,*999)
                  CALL field_parameter_set_update_finish(dependent_field,variable_type,field_analytic_values_set_type, &
                    & err,error,*999)
                ELSE
                  CALL flagerror("Field variable is not associated.",err,error,*999)
                ENDIF
              ENDDO
              CALL field_parameter_set_data_restore(geometric_field,field_u_variable_type,field_values_set_type, &
                & geometric_parameters,err,error,*999)
            ELSE
              CALL flagerror("Boundary conditions is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Equations set geometric field is not associated.",err,error,*999)
          ENDIF            
        ELSE
          CALL flagerror("Equations set dependent field is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Equations set analytic is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
    
    exits("Poiseuille_BoundaryConditionsAnalyticCalculate")
    RETURN
999 errors("Poiseuille_BoundaryConditionsAnalyticCalculate",err,error)
    exits("Poiseuille_BoundaryConditionsAnalyticCalculate")
    RETURN 1
    
  END SUBROUTINE poiseuille_boundaryconditionsanalyticcalculate
  
  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_equation_equations_set_setup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("POISEUILLE_EQUATION_EQUATIONS_SET_SETUP",err,error,*999)

    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a Poiseuille type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_static_poiseuille_subtype)
        CALL poiseuille_equationssetstaticsetup(equations_set,equations_set_setup,err,error,*999)
      CASE(equations_set_dynamic_poiseuille_subtype)
        CALL poiseuille_equationssetstaticsetup(equations_set,equations_set_setup,err,error,*999)
      CASE DEFAULT
        local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Poiseuille equation type of a fluid mechanics equation set class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("POISEUILLE_EQUATION_EQUATIONS_SET_SETUP")
    RETURN
999 errorsexits("POISEUILLE_EQUATION_EQUATIONS_SET_SETUP",err,error)
    RETURN 1
  END SUBROUTINE poiseuille_equation_equations_set_setup

  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_equationssetsolutionmethodset(EQUATIONS_SET,SOLUTION_METHOD,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: SOLUTION_METHOD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("Poiseuille_EquationsSetSolutionMethodSet",err,error,*999)
    
    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a Poiseuille type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_static_poiseuille_subtype,equations_set_dynamic_poiseuille_subtype)        
        SELECT CASE(solution_method)
        CASE(equations_set_fem_solution_method)
          equations_set%SOLUTION_METHOD=equations_set_fem_solution_method
        CASE(equations_set_bem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fd_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_fv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfem_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE(equations_set_gfv_solution_method)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          local_error="The specified solution method of "//trim(number_to_vstring(solution_method,"*",err,error))//" is invalid."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      CASE DEFAULT
        local_error="Equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Poiseuille equation type of an fluid mechanics equations set class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Poiseuille_EquationsSetSolutionMethodSet")
    RETURN
999 errorsexits("Poiseuille_EquationsSetSolutionMethodSet",err,error)
    RETURN 1
    
  END SUBROUTINE poiseuille_equationssetsolutionmethodset

  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_equationssetspecificationset(equationsSet,specification,err,error,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: equationsSet
    INTEGER(INTG), INTENT(IN) :: specification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: localError
    INTEGER(INTG) :: subtype

    enters("Poiseuille_EquationsSetSpecificationSet",err,error,*999)

    IF(ASSOCIATED(equationsset)) THEN
      IF(SIZE(specification,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a Poiseuille equation type equations set.", &
          & err,error,*999)
      END IF
      subtype=specification(3)
      SELECT CASE(subtype)
      CASE(equations_set_static_poiseuille_subtype, &
          & equations_set_dynamic_poiseuille_subtype)
        !ok
      CASE DEFAULT
        localerror="The third equations set specification of "//trim(numbertovstring(subtype,"*",err,error))// &
          & " is not valid for a Poiseuille equation type of a fluid mechanics equations set class."
        CALL flagerror(localerror,err,error,*999)
      END SELECT
      !Set full specification
      IF(ALLOCATED(equationsset%specification)) THEN
        CALL flagerror("Equations set specification is already allocated.",err,error,*999)
      ELSE
        ALLOCATE(equationsset%specification(3),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate equations set specification.",err,error,*999)
      END IF
      equationsset%specification(1:3)=[equations_set_fluid_mechanics_class,equations_set_poiseuille_equation_type,subtype]
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    END IF

    exits("Poiseuille_EquationsSetSpecificationSet")
    RETURN
999 errors("Poiseuille_EquationsSetSpecificationSet",err,error)
    exits("Poiseuille_EquationsSetSpecificationSet")
    RETURN 1
    
  END SUBROUTINE poiseuille_equationssetspecificationset


  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_equationssetstaticsetup(EQUATIONS_SET,EQUATIONS_SET_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    TYPE(equations_set_setup_type), INTENT(INOUT) :: EQUATIONS_SET_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: component_idx,GEOMETRIC_COMPONENT_NUMBER,GEOMETRIC_SCALING_TYPE
    TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("Poiseuille_EquationsSetStaticSetup",err,error,*999)

    NULLIFY(equations)
    NULLIFY(equations_mapping)
    NULLIFY(equations_matrices)
    NULLIFY(geometric_decomposition)
    
    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a Poiseuille type equations set.", &
          & err,error,*999)
      END IF
      IF(equations_set%SPECIFICATION(3)==equations_set_static_poiseuille_subtype) THEN
        SELECT CASE(equations_set_setup%SETUP_TYPE)
        CASE(equations_set_setup_initial_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            CALL poiseuille_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method,err, &
              & error,*999)
          CASE(equations_set_setup_finish_action)
!!TODO: Check valid setup
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_geometry_type)
          !Do nothing???
        CASE(equations_set_setup_dependent_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
              !Create the auto created dependent field
              CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_set%DEPENDENT% &
                & dependent_field,err,error,*999)
              CALL field_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_dependent_type,err,error,*999)
              CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
              CALL field_mesh_decomposition_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_decomposition, &
                & err,error,*999)
              CALL field_geometric_field_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,equations_set%GEOMETRY% &
                & geometric_field,err,error,*999)
              CALL field_number_of_variables_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,2,err,error,*999)
              CALL field_variable_types_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,[field_u_variable_type, &
                & field_deludeln_variable_type],err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_vector_dimension_type,err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_vector_dimension_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_dp_type,err,error,*999)
              CALL field_data_type_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_dp_type,err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,2, &
                & err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,2, &
                & err,error,*999)
              !Default to the geometric interpolation setup
              CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type,1, &
                & geometric_component_number,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                & geometric_component_number,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,2, &
                & geometric_component_number,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & geometric_component_number,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,2, &
                & geometric_component_number,err,error,*999)
              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                  & field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,2, &
                  & field_element_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_deludeln_variable_type,1,field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD, &
                  & field_deludeln_variable_type,2,field_element_based_interpolation,err,error,*999)
                !Default the scaling to the geometric field scaling
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_set%DEPENDENT%DEPENDENT_FIELD,geometric_scaling_type,err,error,*999)
              CASE(equations_set_bem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fd_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ELSE
              !Check the user specified field
              CALL field_type_check(equations_set_setup%FIELD,field_general_type,err,error,*999)
              CALL field_dependent_type_check(equations_set_setup%FIELD,field_dependent_type,err,error,*999)
              CALL field_number_of_variables_check(equations_set_setup%FIELD,2,err,error,*999)
              CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type,field_deludeln_variable_type], &
                & err,error,*999)
              CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type,err,error,*999)
              CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_vector_dimension_type, &
                & err,error,*999)
              CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
              CALL field_data_type_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_dp_type,err,error,*999)
              CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,2,err,error,*999)
              CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,2,err,error,*999)
              SELECT CASE(equations_set%SOLUTION_METHOD)
              CASE(equations_set_fem_solution_method)
                CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,1, &
                  & field_node_based_interpolation,err,error,*999)
                CALL field_component_interpolation_check(equations_set_setup%FIELD,field_u_variable_type,2, &
                  & field_element_based_interpolation,err,error,*999)
                CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,1, &
                  & field_node_based_interpolation,err,error,*999)
                 CALL field_component_interpolation_check(equations_set_setup%FIELD,field_deludeln_variable_type,2, &
                  & field_element_based_interpolation,err,error,*999)
              CASE(equations_set_bem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fd_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_fv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfem_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE(equations_set_gfv_solution_method)
                CALL flagerror("Not implemented.",err,error,*999)
              CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                  & " is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
            ENDIF
          CASE(equations_set_setup_finish_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FIELD_AUTO_CREATED) THEN
              CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
            ENDIF
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation"
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_materials_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            equations_materials=>equations_set%MATERIALS
            IF(ASSOCIATED(equations_materials)) THEN
              IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
                !Create the auto created materials field
                CALL field_create_start(equations_set_setup%FIELD_USER_NUMBER,equations_set%REGION,equations_materials% &
                  & materials_field,err,error,*999)
                CALL field_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_material_type,err,error,*999)
                CALL field_dependent_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_independent_type,err,error,*999)
                CALL field_mesh_decomposition_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_decomposition,err,error,*999)
                CALL field_mesh_decomposition_set_and_lock(equations_materials%MATERIALS_FIELD,geometric_decomposition, &
                  & err,error,*999)
                CALL field_geometric_field_set_and_lock(equations_materials%MATERIALS_FIELD,equations_set%GEOMETRY% &
                  & geometric_field,err,error,*999)
                CALL field_number_of_variables_set_and_lock(equations_materials%MATERIALS_FIELD,1,err,error,*999)
                CALL field_variable_types_set_and_lock(equations_materials%MATERIALS_FIELD,[field_u_variable_type], &
                  & err,error,*999)
                CALL field_dimension_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_vector_dimension_type,err,error,*999)
                CALL field_data_type_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_dp_type,err,error,*999)
                !Set the number of materials components
                CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & 3,err,error,*999)
                !Default the 3 materials components to the first geometric interpolation setup with constant interpolation
                CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & 1,geometric_component_number,err,error,*999)
                DO component_idx=1,3
                  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,geometric_component_number,err,error,*999)
                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & component_idx,field_constant_interpolation,err,error,*999)
                ENDDO !component_idx
                !Default the field scaling to that of the geometric field
                CALL field_scaling_type_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,geometric_scaling_type,err,error,*999)
                CALL field_scaling_type_set(equations_materials%MATERIALS_FIELD,geometric_scaling_type,err,error,*999)
              ELSE
                !Check the user specified field
                CALL field_type_check(equations_set_setup%FIELD,field_material_type,err,error,*999)
                CALL field_dependent_type_check(equations_set_setup%FIELD,field_independent_type,err,error,*999)
                CALL field_number_of_variables_check(equations_set_setup%FIELD,1,err,error,*999)
                CALL field_variable_types_check(equations_set_setup%FIELD,[field_u_variable_type],err,error,*999)
                CALL field_dimension_check(equations_set_setup%FIELD,field_u_variable_type,field_vector_dimension_type, &
                  & err,error,*999)
                CALL field_data_type_check(equations_set_setup%FIELD,field_u_variable_type,field_dp_type,err,error,*999)
                CALL field_number_of_components_check(equations_set_setup%FIELD,field_u_variable_type,3,err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Equations set materials is not associated.",err,error,*999)
            ENDIF

          CASE(equations_set_setup_finish_action)
            equations_materials=>equations_set%MATERIALS
            IF(ASSOCIATED(equations_materials)) THEN
              IF(equations_materials%MATERIALS_FIELD_AUTO_CREATED) THEN
                !Finish creating the materials field
                CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
                !Set the default values for the materials field
                !Default component 1 (viscosity) 
                CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_values_set_type,1,0.005_dp,err,error,*999)
                !Default component 2 (radius) 
                CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_values_set_type,2,0.5_dp,err,error,*999)
                !Default component 3 (length) 
                CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & field_values_set_type,3,1.0_dp,err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Equations set materials is not associated.",err,error,*999)
            ENDIF
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_source_type)
          SELECT CASE(equations_set%SPECIFICATION(3))
          CASE(equations_set_static_poiseuille_subtype)
              SELECT CASE(equations_set_setup%ACTION_TYPE)
               !Set start action
                CASE(equations_set_setup_start_action)
                  !Add in gravity source field
                !Specify finish action
                CASE(equations_set_setup_finish_action)
                  !Add in gravity source field
                CASE DEFAULT
                  local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
                  & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
                  & " is invalid for a linear Poiseuille subtype"
                  CALL flagerror(local_error,err,error,*999)
              END SELECT
            CASE DEFAULT
              local_error="The equation set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
                & " is invalid for a linear Poiseuille equation."
               CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_analytic_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)

          CASE(equations_set_setup_finish_action)
            
          CASE(equations_set_setup_generate_action)
             
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_equations_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) THEN
              !Create the equations
              CALL equations_create_start(equations_set,equations,err,error,*999)
              CALL equations_linearity_type_set(equations,equations_linear,err,error,*999)
              CALL equations_time_dependence_type_set(equations,equations_static,err,error,*999)
            ELSE
              CALL flagerror("Equations set dependent field has not been finished.",err,error,*999)
            ENDIF
          CASE(equations_set_setup_finish_action)
            SELECT CASE(equations_set%SOLUTION_METHOD)
            CASE(equations_set_fem_solution_method)
              !Finish the creation of the equations
              CALL equations_set_equations_get(equations_set,equations,err,error,*999)
              CALL equations_create_finish(equations,err,error,*999)
              !Create the equations mapping.
              CALL equations_mapping_create_start(equations,equations_mapping,err,error,*999)
              CALL equationsmapping_linearmatricesnumberset(equations_mapping,1,err,error,*999)
              CALL equationsmapping_linearmatricesvariabletypesset(equations_mapping,[field_u_variable_type], &
                & err,error,*999)
              CALL equations_mapping_rhs_variable_type_set(equations_mapping,field_deludeln_variable_type,err,error,*999)
              CALL equations_mapping_create_finish(equations_mapping,err,error,*999)
              !Create the equations matrices
              CALL equations_matrices_create_start(equations,equations_matrices,err,error,*999)
              SELECT CASE(equations%SPARSITY_TYPE)
              CASE(equations_matrices_full_matrices)
                CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_block_storage_type], &
                  & err,error,*999)
              CASE(equations_matrices_sparse_matrices)
                CALL equations_matrices_linear_storage_type_set(equations_matrices,[matrix_compressed_row_storage_type], &
                  & err,error,*999)
                CALL equationsmatrices_linearstructuretypeset(equations_matrices,[equations_matrix_fem_structure], &
                  & err,error,*999)
              CASE DEFAULT
                local_error="The equations matrices sparsity type of "// &
                  & trim(number_to_vstring(equations%SPARSITY_TYPE,"*",err,error))//" is invalid."
                CALL flagerror(local_error,err,error,*999)
              END SELECT
              CALL equations_matrices_create_finish(equations_matrices,err,error,*999)
            CASE(equations_set_bem_solution_method)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(equations_set_fd_solution_method)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(equations_set_fv_solution_method)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(equations_set_gfem_solution_method)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE(equations_set_gfv_solution_method)
              CALL flagerror("Not implemented.",err,error,*999)
            CASE DEFAULT
                local_error="The solution method of "//trim(number_to_vstring(equations_set%SOLUTION_METHOD,"*",err,error))// &
                & " is invalid."
              CALL flagerror(local_error,err,error,*999)
            END SELECT
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(equations_set_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The setup type of "//trim(number_to_vstring(equations_set_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a static Poiseuille equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        local_error="The equations set subtype of "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
          & " is not a static Poiseuille equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Poiseuille_EquationsSetStaticSetup")
    RETURN
999 errorsexits("Poiseuille_EquationsSetStaticSetup",err,error)
    RETURN 1
    
  END SUBROUTINE poiseuille_equationssetstaticsetup

  !
  !================================================================================================================================
  !
 
  SUBROUTINE poiseuille_equation_problem_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("POISEUILLE_EQUATION_PROBLEM_SETUP",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(problem%specification)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(problem%specification,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a Poiseuille problem.",err,error,*999)
      END IF
      SELECT CASE(problem%SPECIFICATION(3))
      CASE(problem_static_poiseuille_subtype)
        CALL poiseuille_equation_problem_static_setup(problem,problem_setup,err,error,*999)
      CASE(problem_dynamic_poiseuille_subtype)
        CALL poiseuille_equation_problem_static_setup(problem,problem_setup,err,error,*999)
      CASE DEFAULT
        local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Poiseuille equation type of a fluid mechanics problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("POISEUILLE_EQUATION_PROBLEM_SETUP")
    RETURN
999 errorsexits("POISEUILLE_EQUATION_PROBLEM_SETUP",err,error)
    RETURN 1
  END SUBROUTINE poiseuille_equation_problem_setup
  
  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_finiteelementcalculate(EQUATIONS_SET,ELEMENT_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(equations_set_type), POINTER :: EQUATIONS_SET
    INTEGER(INTG), INTENT(IN) :: ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) FIELD_VAR_TYPE,ng
    TYPE(basis_type), POINTER :: DEPENDENT_BASIS,GEOMETRIC_BASIS
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_mapping_linear_type), POINTER :: LINEAR_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_matrices_linear_type), POINTER :: LINEAR_MATRICES
    TYPE(equations_matrices_rhs_type), POINTER :: RHS_VECTOR
    TYPE(equations_matrix_type), POINTER :: EQUATIONS_MATRIX
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE
    TYPE(quadrature_scheme_type), POINTER :: QUADRATURE_SCHEME
    TYPE(varying_string) :: LOCAL_ERROR

    enters("Poiseuille_FiniteElementCalculate",err,error,*999)

    IF(ASSOCIATED(equations_set)) THEN
      IF(.NOT.ALLOCATED(equations_set%SPECIFICATION)) THEN
        CALL flagerror("Equations set specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(equations_set%SPECIFICATION,1)/=3) THEN
        CALL flagerror("Equations set specification must have three entries for a Poiseuille type equations set.", &
          & err,error,*999)
      END IF
      equations=>equations_set%EQUATIONS
      IF(ASSOCIATED(equations)) THEN
        SELECT CASE(equations_set%SPECIFICATION(3))
        CASE(equations_set_static_poiseuille_subtype,equations_set_dynamic_poiseuille_subtype)
          !Store all these in equations matrices/somewhere else?????
          dependent_field=>equations%INTERPOLATION%DEPENDENT_FIELD
          geometric_field=>equations%INTERPOLATION%GEOMETRIC_FIELD
          materials_field=>equations%INTERPOLATION%MATERIALS_FIELD
          equations_matrices=>equations%EQUATIONS_MATRICES
          linear_matrices=>equations_matrices%LINEAR_MATRICES
          equations_matrix=>linear_matrices%MATRICES(1)%PTR
          rhs_vector=>equations_matrices%RHS_VECTOR
          equations_mapping=>equations%EQUATIONS_MAPPING
          linear_mapping=>equations_mapping%LINEAR_MAPPING
          field_variable=>linear_mapping%EQUATIONS_MATRIX_TO_VAR_MAPS(1)%VARIABLE
          field_var_type=field_variable%VARIABLE_TYPE
          dependent_basis=>dependent_field%DECOMPOSITION%DOMAIN(dependent_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
            & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
          geometric_basis=>geometric_field%DECOMPOSITION%DOMAIN(geometric_field%DECOMPOSITION%MESH_COMPONENT_NUMBER)%PTR% &
            & topology%ELEMENTS%ELEMENTS(element_number)%BASIS
          quadrature_scheme=>dependent_basis%QUADRATURE%QUADRATURE_SCHEME_MAP(basis_default_quadrature_scheme)%PTR
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
            & geometric_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,equations%INTERPOLATION% &
            & materials_interp_parameters(field_u_variable_type)%PTR,err,error,*999)
          !Loop over gauss points
          DO ng=1,quadrature_scheme%NUMBER_OF_GAUSS
            CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & geometric_interp_point(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolate_gauss(no_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
            CALL field_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
              & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)

!TODO

            
          ENDDO !ng
        CASE DEFAULT
          local_error="Equations set subtype "//trim(number_to_vstring(equations_set%SPECIFICATION(3),"*",err,error))// &
            & " is not valid for a Poiseuille equation type of a fluid mechanics equations set class."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        CALL flagerror("Equations set equations is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Poiseuille_FiniteElementCalculate")
    RETURN
999 errorsexits("Poiseuille_FiniteElementCalculate",err,error)
    RETURN 1
    
  END SUBROUTINE poiseuille_finiteelementcalculate

  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_problemspecificationset(problem,problemSpecification,err,error,*)

    !Argument variables
    TYPE(problem_type), POINTER :: problem
    INTEGER(INTG), INTENT(IN) :: problemSpecification(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: localError
    INTEGER(INTG) :: problemSubtype

    enters("Poiseuille_ProblemSpecificationSet",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(SIZE(problemspecification,1)==3) THEN
        problemsubtype=problemspecification(3)
        SELECT CASE(problemsubtype)
        CASE(problem_static_poiseuille_subtype, &
            & problem_dynamic_poiseuille_subtype)
          !ALl ok
        CASE DEFAULT
          localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
            & " is not valid for a Poiseuille fluid mechanics problem."
          CALL flagerror(localerror,err,error,*999)
        END SELECT
        IF(ALLOCATED(problem%specification)) THEN
          CALL flagerror("Problem specification is already allocated.",err,error,*999)
        ELSE
          ALLOCATE(problem%specification(3),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate problem specification.",err,error,*999)
        END IF
        problem%specification(1:3)=[problem_fluid_mechanics_class,problem_poiseuille_equation_type,problemsubtype]
      ELSE
        CALL flagerror("Poiseuille problem specification must have three entries.",err,error,*999)
      END IF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    END IF

    exits("Poiseuille_ProblemSpecificationSet")
    RETURN
999 errors("Poiseuille_ProblemSpecificationSet",err,error)
    exits("Poiseuille_ProblemSpecificationSet")
    RETURN 1
    
  END SUBROUTINE poiseuille_problemspecificationset

  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_equation_problem_static_setup(PROBLEM,PROBLEM_SETUP,ERR,ERROR,*)

    !Argument variables
    TYPE(problem_type), POINTER :: PROBLEM
    TYPE(problem_setup_type), INTENT(INOUT) :: PROBLEM_SETUP
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP,CONTROL_LOOP_ROOT
    TYPE(solver_type), POINTER :: SOLVER
    TYPE(solver_equations_type), POINTER :: SOLVER_EQUATIONS
    TYPE(solvers_type), POINTER :: SOLVERS
    TYPE(varying_string) :: LOCAL_ERROR
   
    enters("POISEUILLE_EQUATION_PROBLEM_STATIC_SETUP",err,error,*999)

    NULLIFY(control_loop)
    NULLIFY(solver)
    NULLIFY(solver_equations)
    NULLIFY(solvers)
    IF(ASSOCIATED(problem)) THEN
      IF(.NOT.ALLOCATED(problem%specification)) THEN
        CALL flagerror("Problem specification is not allocated.",err,error,*999)
      ELSE IF(SIZE(problem%specification,1)<3) THEN
        CALL flagerror("Problem specification must have three entries for a Poiseuille problem.",err,error,*999)
      END IF
      IF(problem%SPECIFICATION(3)==problem_static_poiseuille_subtype) THEN
        SELECT CASE(problem_setup%SETUP_TYPE)
        CASE(problem_setup_initial_type)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Do nothing????
          CASE(problem_setup_finish_action)
            !Do nothing????
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_control_type)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Set up a simple control loop
            CALL control_loop_create_start(problem,control_loop,err,error,*999)
          CASE(problem_setup_finish_action)
            !Finish the control loops
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            CALL control_loop_create_finish(control_loop,err,error,*999)            
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_solvers_type)
          !Get the control loop
          control_loop_root=>problem%CONTROL_LOOP
          CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Start the solvers creation
            CALL solvers_create_start(control_loop,solvers,err,error,*999)
            CALL solvers_number_set(solvers,1,err,error,*999)
            !Set the solver to be a linear solver
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
             !Start the linear solver creation
            CALL solver_type_set(solver,solver_linear_type,err,error,*999)
            !Set solver defaults
            CALL solver_library_type_set(solver,solver_petsc_library,err,error,*999)
          CASE(problem_setup_finish_action)
            !Get the solvers
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            !Finish the solvers creation
            CALL solvers_create_finish(solvers,err,error,*999)
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(problem_setup_solver_equations_type)
          SELECT CASE(problem_setup%ACTION_TYPE)
          CASE(problem_setup_start_action)
            !Get the control loop
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            !Get the solver
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            !Create the solver equations
            CALL solver_equations_create_start(solver,solver_equations,err,error,*999)
            CALL solver_equations_linearity_type_set(solver_equations,solver_equations_linear,err,error,*999)
            CALL solver_equations_time_dependence_type_set(solver_equations,solver_equations_static,err,error,*999)
            CALL solver_equations_sparsity_type_set(solver_equations,solver_sparse_matrices,err,error,*999)
          CASE(problem_setup_finish_action)
            !Get the control loop
            control_loop_root=>problem%CONTROL_LOOP
            CALL control_loop_get(control_loop_root,control_loop_node,control_loop,err,error,*999)
            !Get the solver equations
            CALL control_loop_solvers_get(control_loop,solvers,err,error,*999)
            CALL solvers_solver_get(solvers,1,solver,err,error,*999)
            CALL solver_solver_equations_get(solver,solver_equations,err,error,*999)
            !Finish the solver equations creation
            CALL solver_equations_create_finish(solver_equations,err,error,*999)             
          CASE DEFAULT
            local_error="The action type of "//trim(number_to_vstring(problem_setup%ACTION_TYPE,"*",err,error))// &
              & " for a setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
              & " is invalid for a static Poiseuille equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE DEFAULT
          local_error="The setup type of "//trim(number_to_vstring(problem_setup%SETUP_TYPE,"*",err,error))// &
            & " is invalid for a static Poiseuille equation."
          CALL flagerror(local_error,err,error,*999)
        END SELECT
      ELSE
        local_error="The problem subtype of "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
          & " does not equal a static Poiseuille equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("POISEUILLE_EQUATION_PROBLEM_STATIC_SETUP")
    RETURN
999 errorsexits("POISEUILLE_EQUATION_PROBLEM_STATIC_SETUP",err,error)
    RETURN 1
  END SUBROUTINE poiseuille_equation_problem_static_setup
  
  !
  !================================================================================================================================
  !

  SUBROUTINE poiseuille_post_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(solver_type), POINTER :: SOLVER2
    TYPE(varying_string) :: LOCAL_ERROR

    enters("POISEUILLE_POST_SOLVE",err,error,*999)
    NULLIFY(solver2)
    IF(ASSOCIATED(control_loop)) THEN
      IF(ASSOCIATED(solver)) THEN
        IF(ASSOCIATED(control_loop%PROBLEM)) THEN 
          IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
            CALL flagerror("Problem specification is not allocated.",err,error,*999)
          ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
            CALL flagerror("Problem specification must have three entries for a Poiseuille problem.",err,error,*999)
          END IF
          SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
            CASE(problem_static_poiseuille_subtype)
!               do nothing
            CASE(problem_dynamic_poiseuille_subtype)
!               do nothing
            CASE DEFAULT
              local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
                & " is not valid for a Poiseuille type of a fluid mechanics problem class."
              CALL flagerror(local_error,err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("Problem is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Solver is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Control loop is not associated.",err,error,*999)
    ENDIF

    exits("POISEUILLE_POST_SOLVE")
    RETURN
999 errorsexits("POISEUILLE_POST_SOLVE",err,error)
    RETURN 1
  END SUBROUTINE poiseuille_post_solve

  !
  !================================================================================================================================
  !


  SUBROUTINE poiseuille_pre_solve(CONTROL_LOOP,SOLVER,ERR,ERROR,*)

    !Argument variables
    TYPE(control_loop_type), POINTER :: CONTROL_LOOP
    TYPE(solver_type), POINTER :: SOLVER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(solver_type), POINTER :: SOLVER2
    TYPE(varying_string) :: LOCAL_ERROR

    enters("POISEUILLE_PRE_SOLVE",err,error,*999)
    NULLIFY(solver2)
 
    IF(ASSOCIATED(control_loop)) THEN
      IF(ASSOCIATED(solver)) THEN
        IF(ASSOCIATED(control_loop%PROBLEM)) THEN
          IF(.NOT.ALLOCATED(control_loop%problem%specification)) THEN
            CALL flagerror("Problem specification is not allocated.",err,error,*999)
          ELSE IF(SIZE(control_loop%problem%specification,1)<3) THEN
            CALL flagerror("Problem specification must have three entries for a Poiseuille problem.",err,error,*999)
          END IF
          SELECT CASE(control_loop%PROBLEM%SPECIFICATION(3))
            CASE(problem_static_poiseuille_subtype)
!               do nothing
            CASE(problem_dynamic_poiseuille_subtype)
!               do nothing
            CASE DEFAULT
              local_error="Problem subtype "//trim(number_to_vstring(control_loop%PROBLEM%SPECIFICATION(3),"*",err,error))// &
                & " is not valid for a Poiseuille type of a fluid mechanics problem class."
              CALL flagerror(local_error,err,error,*999)
          END SELECT
        ELSE
          CALL flagerror("Problem is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Solver is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Control loop is not associated.",err,error,*999)
    ENDIF

    exits("POISEUILLE_PRE_SOLVE")
    RETURN
999 errorsexits("POISEUILLE_PRE_SOLVE",err,error)
    RETURN 1
  END SUBROUTINE poiseuille_pre_solve

  !
  !================================================================================================================================
  !
 

END MODULE poiseuille_equations_routines