Helmholtz_equations_routines.f90

Go to the documentation of this file.

MODULE helmholtz_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 node_routines
  USE problem_constants
  USE strings
  USE solver_routines
  USE timer
  USE types

#include "macros.h"  

  IMPLICIT NONE

  PRIVATE

  !Module parameters

  !Module types

  !Module variables

  !Interfaces

!!MERGE: move

  PUBLIC helmholtz_boundaryconditionsanalyticcalculate
  
  PUBLIC helmholtz_equationssetsolutionmethodset
  
  PUBLIC helmholtz_equation_equations_set_setup

  PUBLIC helmholtz_equationssetspecificationset
  
  PUBLIC helmholtz_equation_finite_element_calculate

  PUBLIC helmholtz_equation_problem_setup

  PUBLIC helmholtz_problemspecificationset

CONTAINS

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


  SUBROUTINE helmholtz_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,node_idx,NUMBER_OF_DIMENSIONS,variable_idx,variable_type,local_ny,deriv_idx,dim_idx
    REAL(DP), POINTER :: GEOMETRIC_PARAMETERS(:), MATERIALS_PARAMETERS(:)
    TYPE(domain_type), POINTER :: DOMAIN
    TYPE(domain_nodes_type), POINTER :: DOMAIN_NODES
    TYPE(field_type), POINTER :: DEPENDENT_FIELD,GEOMETRIC_FIELD,MATERIALS_FIELD
    TYPE(field_variable_type), POINTER :: FIELD_VARIABLE,GEOMETRIC_VARIABLE,MATERIALS_VARIABLE
    REAL(DP) :: VALUE,X(3),k,mu
    TYPE(varying_string) :: LOCAL_ERROR    
    
    enters("Helmholtz_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            
            materials_field=>equations_set%MATERIALS%MATERIALS_FIELD
            IF(ASSOCIATED(materials_field)) THEN            
              CALL field_number_of_components_get(geometric_field,field_u_variable_type,number_of_dimensions,err,error,*999)
              NULLIFY(geometric_variable)
              NULLIFY(geometric_parameters)
              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)
              NULLIFY(materials_variable)
              NULLIFY(materials_parameters)
              CALL field_variable_get(materials_field,field_u_variable_type,materials_variable,err,error,*999)
              CALL field_parameter_set_data_get(materials_field,field_u_variable_type,field_values_set_type,materials_parameters, &
                & err,error,*999)
              IF(ASSOCIATED(boundary_conditions)) THEN
                DO variable_idx=1,dependent_field%NUMBER_OF_VARIABLES
                  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
                      IF(field_variable%COMPONENTS(component_idx)%INTERPOLATION_TYPE==field_node_based_interpolation) THEN
                        IF(materials_variable%COMPONENTS(1)%INTERPOLATION_TYPE==field_constant_interpolation) THEN
                          local_ny=materials_variable%COMPONENTS(1)%PARAM_TO_DOF_MAP%CONSTANT_PARAM2DOF_MAP
                          k=materials_parameters(local_ny)
                          mu=k/sqrt(2.0_dp)
                          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
                                    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_helmholtz_equation_two_dim_1)
                                      !u=cos(k*x/sqrt(2))*sin(k*y/sqrt(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=cos(mu*x(1))*sin(mu*x(2))
                                        CASE(global_deriv_s1)
                                          VALUE=-mu*sin(mu*x(1))*sin(mu*x(2))
                                        CASE(global_deriv_s2)
                                          VALUE=mu*cos(mu*x(1))*cos(mu*x(2))
                                        CASE(global_deriv_s1_s2)
                                          VALUE=-mu*mu*sin(mu*x(1))*cos(mu*x(2))
                                        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)
                                          VALUE=0.0_dp !!TODO
                                        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 DEFAULT
                                        local_error="The variable type of "//trim(number_to_vstring(variable_type,"*",err,error))//&
                                          & " is invalid."
                                        CALL flagerror(local_error,err,error,*999)
                                      END SELECT
                                    CASE DEFAULT
                                      local_error="The analytic function type of "// &
                                        & trim(number_to_vstring(equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
                                        & " is invalid."
                                      CALL flagerror(local_error,err,error,*999)
                                    END SELECT
                                    local_ny=field_variable%COMPONENTS(component_idx)%PARAM_TO_DOF_MAP% &
                                      & node_param2dof_map%NODES(node_idx)%DERIVATIVES(deriv_idx)%VERSIONS(1)
                                    CALL field_parameter_set_update_local_dof(dependent_field,variable_type, &
                                      & field_analytic_values_set_type,local_ny,VALUE,err,error,*999)
                                    IF(variable_type==field_u_variable_type) 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,dependent_field,variable_type, &
                                          & local_ny,boundary_condition_fixed,VALUE,err,error,*999)
                                      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 constant interpolation of the material field is implemented.",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 !variable_idx
                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 materials field 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("Helmholtz_BoundaryConditionsAnalyticCalculate")
    RETURN
999 errors("Helmholtz_BoundaryConditionsAnalyticCalculate",err,error)
    exits("Helmholtz_BoundaryConditionsAnalyticCalculate")
    RETURN 1
    
  END SUBROUTINE helmholtz_boundaryconditionsanalyticcalculate
  
  !
  !================================================================================================================================
  !

  SUBROUTINE helmholtz_equation_finite_element_calculate(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,mh,mhs,mi,ms,nh,nhs,ni,ns
    REAL(DP) :: RWG,SUM,PGMSI(3),PGNSI(3),PGM,PGN,k
    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
    
#ifdef TAUPROF
    CHARACTER(26) :: CVAR
    INTEGER :: GAUSS_POINT_LOOP_PHASE(2) = (/ 0, 0 /)
    SAVE gauss_point_loop_phase
#endif

    enters("HELMHOLTZ_EQUATION_FINITE_ELEMENT_CALCULATE",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 Helmholtz 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_standard_helmholtz_subtype)
!!TODO: move these and scale factor adjustment out once generalised Helmholtz is put in.
          !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
#ifdef TAUPROF
              WRITE (cvar,'(a17,i2)') 'Gauss Point Loop ',ng
              CALL tau_phase_create_dynamic(gauss_point_loop_phase,cvar)
              CALL tau_phase_start(gauss_point_loop_phase)
#endif
            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_interpolated_point_metrics_calculate(geometric_basis%NUMBER_OF_XI,equations%INTERPOLATION% &
              & geometric_interp_point_metrics(field_u_variable_type)%PTR,err,error,*999)

            CALL field_interpolate_gauss(first_part_deriv,basis_default_quadrature_scheme,ng,equations%INTERPOLATION% &
              & materials_interp_point(field_u_variable_type)%PTR,err,error,*999)
            k = equations%INTERPOLATION%MATERIALS_INTERP_POINT(field_u_variable_type)%PTR%VALUES(1,no_part_deriv) ! wave number

            !Calculate RWG.
!!TODO: Think about symmetric problems. 
            rwg=equations%INTERPOLATION%GEOMETRIC_INTERP_POINT_METRICS(field_u_variable_type)%PTR%JACOBIAN* &
              & quadrature_scheme%GAUSS_WEIGHTS(ng)
            !Loop over field components
            mhs=0          
            DO mh=1,field_variable%NUMBER_OF_COMPONENTS
              !Loop over element rows
!!TODO: CHANGE ELEMENT CALCULATE TO WORK OF ns ???
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                mhs=mhs+1
                nhs=0
                IF(equations_matrix%UPDATE_MATRIX) THEN
                  !Loop over element columns
                  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                    DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                      nhs=nhs+1
                      pgm=quadrature_scheme%GAUSS_BASIS_FNS(ms,no_part_deriv,ng)
                      pgn=quadrature_scheme%GAUSS_BASIS_FNS(ns,no_part_deriv,ng)
                      DO ni=1,dependent_basis%NUMBER_OF_XI
                        pgmsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ms,partial_derivative_first_derivative_map(ni),ng)
                        pgnsi(ni)=quadrature_scheme%GAUSS_BASIS_FNS(ns,partial_derivative_first_derivative_map(ni),ng)
                      ENDDO !ni
                      sum=0.0_dp
                      sum=sum-k*k*pgm*pgn
                      DO mi=1,dependent_basis%NUMBER_OF_XI
                        DO ni=1,dependent_basis%NUMBER_OF_XI
                          sum=sum+pgmsi(mi)*pgnsi(ni)*equations%INTERPOLATION% &
                            & geometric_interp_point_metrics(field_u_variable_type)%PTR%GU(mi,ni)
                        ENDDO !ni
                      ENDDO !mi
                      equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)+sum*rwg
                    ENDDO !ns
                  ENDDO !nh
                ENDIF
                IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=0.0_dp
              ENDDO !ms
            ENDDO !mh
#ifdef TAUPROF
            CALL tau_phase_stop(gauss_point_loop_phase)
#endif
          ENDDO !ng
          
          !Scale factor adjustment
          IF(dependent_field%SCALINGS%SCALING_TYPE/=field_no_scaling) THEN
            CALL field_interpolationparametersscalefactorselementget(element_number,equations%INTERPOLATION% &
              & dependent_interp_parameters(field_var_type)%PTR,err,error,*999)
            mhs=0          
            DO mh=1,field_variable%NUMBER_OF_COMPONENTS
              !Loop over element rows
              DO ms=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                mhs=mhs+1                    
                nhs=0
                IF(equations_matrix%UPDATE_MATRIX) THEN
                  !Loop over element columns
                  DO nh=1,field_variable%NUMBER_OF_COMPONENTS
                    DO ns=1,dependent_basis%NUMBER_OF_ELEMENT_PARAMETERS
                      nhs=nhs+1
                      equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)=equations_matrix%ELEMENT_MATRIX%MATRIX(mhs,nhs)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)* &
                        & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ns,nh)
                    ENDDO !ns
                  ENDDO !nh
                ENDIF
                IF(rhs_vector%UPDATE_VECTOR) rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)=rhs_vector%ELEMENT_VECTOR%VECTOR(mhs)* &
                  & equations%INTERPOLATION%DEPENDENT_INTERP_PARAMETERS(field_var_type)%PTR%SCALE_FACTORS(ms,mh)
              ENDDO !ms
            ENDDO !mh
          ENDIF
        CASE(equations_set_generalised_helmholtz_subtype)
          CALL flagerror("Not implemented.",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 Helmholtz equation type of a classical field 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("HELMHOLTZ_EQUATION_FINITE_ELEMENT_CALCULATE")
    RETURN
999 errorsexits("HELMHOLTZ_EQUATION_FINITE_ELEMENT_CALCULATE",err,error)
    RETURN 1
  END SUBROUTINE helmholtz_equation_finite_element_calculate

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

  SUBROUTINE helmholtz_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("HELMHOLTZ_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 Helmholtz type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_standard_helmholtz_subtype)
        CALL helmholtz_equationssetstandardsetup(equations_set,equations_set_setup,err,error,*999)
      CASE(equations_set_generalised_helmholtz_subtype)
        CALL flagerror("Not implemented.",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 Helmholtz equation type of a classical field 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("HELMHOLTZ_EQUATION_EQUATIONS_SET_SETUP")
    RETURN
999 errorsexits("HELMHOLTZ_EQUATION_EQUATIONS_SET_SETUP",err,error)
    RETURN 1
  END SUBROUTINE helmholtz_equation_equations_set_setup

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

  SUBROUTINE helmholtz_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("Helmholtz_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 Helmholtz type equations set.", &
          & err,error,*999)
      END IF
      SELECT CASE(equations_set%SPECIFICATION(3))
      CASE(equations_set_standard_helmholtz_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(equations_set_generalised_helmholtz_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 Helmholtz equation type of an classical field 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("Helmholtz_EquationsSetSolutionMethodSet")
    RETURN
999 errors("Helmholtz_EquationsSetSolutionMethodSet",err,error)
    exits("Helmholtz_EquationsSetSolutionMethodSet")
    RETURN 1
    
  END SUBROUTINE helmholtz_equationssetsolutionmethodset

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

  SUBROUTINE helmholtz_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("Helmholtz_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 Helmholtz type equations set.", &
          & err,error,*999)
      END IF
      subtype=specification(3)
      SELECT CASE(subtype)
      CASE(equations_set_standard_helmholtz_subtype)
        !ok
      CASE(equations_set_generalised_helmholtz_subtype)
        CALL flagerror("Not implemented.",err,error,*999)
      CASE DEFAULT
        localerror="The third equations set specification of "//trim(numbertovstring(subtype,"*",err,error))// &
          & " is not valid for a Helmholtz type of a classical field equations set."
        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_classical_field_class,equations_set_helmholtz_equation_type,subtype]
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    END IF

    exits("Helmholtz_EquationsSetSpecificationSet")
    RETURN
999 errors("Helmholtz_EquationsSetSpecificationSet",err,error)
    exits("Helmholtz_EquationsSetSpecificationSet")
    RETURN 1
    
  END SUBROUTINE helmholtz_equationssetspecificationset

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

  SUBROUTINE helmholtz_equationssetstandardsetup(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) :: GEOMETRIC_MESH_COMPONENT,GEOMETRIC_SCALING_TYPE,NUMBER_OF_DIMENSIONS
    TYPE(decomposition_type), POINTER :: GEOMETRIC_DECOMPOSITION
    TYPE(field_type), POINTER :: ANALYTIC_FIELD,DEPENDENT_FIELD,GEOMETRIC_FIELD
    TYPE(equations_type), POINTER :: EQUATIONS
    TYPE(equations_mapping_type), POINTER :: EQUATIONS_MAPPING
    TYPE(equations_matrices_type), POINTER :: EQUATIONS_MATRICES
    TYPE(equations_set_materials_type), POINTER :: EQUATIONS_MATERIALS
    INTEGER(INTG) :: MATERIAL_FIELD_NUMBER_OF_VARIABLES, MATERIAL_FIELD_NUMBER_OF_COMPONENTS
    INTEGER(INTG) :: i,GEOMETRIC_COMPONENT_NUMBER
    TYPE(varying_string) :: LOCAL_ERROR
    
    enters("HELMHOLTZ_EQUATION_EQUATION_SET_STANDARD_SETUP",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 Helmholtz type equations set.", &
          & err,error,*999)
      END IF
      IF(equations_set%SPECIFICATION(3)==equations_set_standard_helmholtz_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 helmholtz_equationssetsolutionmethodset(equations_set,equations_set_fem_solution_method,err,error, &
              & *999)               
          CASE(equations_set_setup_finish_action)
            !Do nothing
          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 standard Helmholtz 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_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,"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_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,"Phi",err,error,*999)
              CALL field_variable_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,"del Phi/del n", &
                & err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type, &
                & field_scalar_dimension_type,err,error,*999)
              CALL field_dimension_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type, &
                & field_scalar_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,1, &
                & err,error,*999)
              CALL field_number_of_components_set_and_lock(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & err,error,*999)
              CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1,"Phi",err,error,*999)
              CALL field_component_label_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & "del Phi/del n",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_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_u_variable_type,1, &
                & geometric_mesh_component,err,error,*999)
              CALL field_component_mesh_component_set(equations_set%DEPENDENT%DEPENDENT_FIELD,field_deludeln_variable_type,1, &
                & geometric_mesh_component,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_deludeln_variable_type,1,field_node_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_scalar_dimension_type,err,error,*999)
              CALL field_dimension_check(equations_set_setup%FIELD,field_deludeln_variable_type,field_scalar_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,1,err,error,*999)
              CALL field_number_of_components_check(equations_set_setup%FIELD,field_deludeln_variable_type,1,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_deludeln_variable_type,1, &
                  & field_node_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 standard Helmholtz equation"
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_materials_type)
          material_field_number_of_variables=1
          material_field_number_of_components=1 ! wave number k
          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(equations_materials%MATERIALS_FIELD, &
                  & material_field_number_of_variables,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)
                CALL field_number_of_components_set_and_lock(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                  & material_field_number_of_components,err,error,*999)
                CALL field_component_mesh_component_get(equations_set%GEOMETRY%GEOMETRIC_FIELD,field_u_variable_type, &
                  & 1,geometric_component_number,err,error,*999)                

                !Default to constant interpolation of material field
                DO i=1,material_field_number_of_components
                  CALL field_component_mesh_component_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & i,geometric_component_number,err,error,*999)
                  CALL field_component_interpolation_set(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & i,field_constant_interpolation,err,error,*999)
                END DO

                !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,1,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
                CALL field_create_finish(equations_materials%MATERIALS_FIELD,err,error,*999)
                !Set the default values for the materials field
                DO i=1,material_field_number_of_components
                  CALL field_component_values_initialise(equations_materials%MATERIALS_FIELD,field_u_variable_type, &
                    & field_values_set_type, i, 1.0_dp, err, error, *999)
                ENDDO
              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 standard Helmholtz equation."
            CALL flagerror(local_error,err,error,*999)
          END SELECT
        CASE(equations_set_setup_source_type)
          SELECT CASE(equations_set_setup%ACTION_TYPE)
          CASE(equations_set_setup_start_action)
            !Do nothing
          CASE(equations_set_setup_finish_action)
            !Do nothing
          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 standard Helmholtz 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)
            IF(equations_set%DEPENDENT%DEPENDENT_FINISHED) 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)
                  SELECT CASE(equations_set_setup%ANALYTIC_FUNCTION_TYPE)
                  CASE(equations_set_helmholtz_equation_two_dim_1)
                    !Check that we are in 2D
                    IF(number_of_dimensions/=2) THEN
                      local_error="The number of geometric dimensions of "// &
                        & trim(number_to_vstring(number_of_dimensions,"*",err,error))// &
                        & " is invalid. The analytic function type of "// &
                        & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
                        & " requires that there be 2 geometric dimensions."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF
                    !Create analytic field if required
                    !Set analtyic function type
                    equations_set%ANALYTIC%ANALYTIC_FUNCTION_TYPE=equations_set_helmholtz_equation_two_dim_1
                  CASE DEFAULT
                    local_error="The specified analytic function type of "// &
                      & trim(number_to_vstring(equations_set_setup%ANALYTIC_FUNCTION_TYPE,"*",err,error))// &
                      & " is invalid for a standard Helmholtz equation."
                    CALL flagerror(local_error,err,error,*999)
                  END SELECT
                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 dependent field has not been finished.",err,error,*999)
            ENDIF
          CASE(equations_set_setup_finish_action)
            IF(ASSOCIATED(equations_set%ANALYTIC)) THEN
              analytic_field=>equations_set%ANALYTIC%ANALYTIC_FIELD
              IF(ASSOCIATED(analytic_field)) THEN
                IF(equations_set%ANALYTIC%ANALYTIC_FIELD_AUTO_CREATED) THEN
                  CALL field_create_finish(equations_set%DEPENDENT%DEPENDENT_FIELD,err,error,*999)
                ENDIF
              ENDIF
            ELSE
              CALL flagerror("Equations set analytic 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 standard Helmholtz 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
              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 equations creation
              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 standard Helmholtz 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 standard Helmholtz 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))// &
          & " does not equal a standard Helmholtz equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Equations set is not associated.",err,error,*999)
    ENDIF
       
    exits("Helmholtz_EquationsSetStandardSetup")
    RETURN
999 errorsexits("Helmholtz_EquationsSetStandardSetup",err,error)
    RETURN 1
    
  END SUBROUTINE helmholtz_equationssetstandardsetup

  !
  !================================================================================================================================
  !
 
  SUBROUTINE helmholtz_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("HELMHOLTZ_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 Helmholtz equation problem.",err,error,*999)
      END IF
      SELECT CASE(problem%SPECIFICATION(3))
      CASE(problem_standard_helmholtz_subtype)
        CALL helmholtz_equation_problem_standard_setup(problem,problem_setup,err,error,*999)
      CASE(problem_generalised_helmholtz_subtype)
        CALL flagerror("Not implemented.",err,error,*999)
      CASE DEFAULT
        local_error="Problem subtype "//trim(number_to_vstring(problem%SPECIFICATION(3),"*",err,error))// &
          & " is not valid for a Helmholtz equation type of a classical field problem class."
        CALL flagerror(local_error,err,error,*999)
      END SELECT
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("HELMHOLTZ_EQUATION_PROBLEM_SETUP")
    RETURN
999 errorsexits("HELMHOLTZ_EQUATION_PROBLEM_SETUP",err,error)
    RETURN 1
  END SUBROUTINE helmholtz_equation_problem_setup
  
  !
  !================================================================================================================================
  !

  SUBROUTINE helmholtz_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

    CALL enters("Helmholtz_ProblemSpecificationSet",err,error,*999)

    IF(ASSOCIATED(problem)) THEN
      IF(SIZE(problemspecification,1)==3) THEN
        problemsubtype=problemspecification(3)
        SELECT CASE(problemsubtype)
        CASE(problem_standard_helmholtz_subtype)
          !ok
        CASE(problem_generalised_helmholtz_subtype)
          CALL flagerror("Not implemented.",err,error,*999)
        CASE DEFAULT
          localerror="The third problem specification of "//trim(numbertovstring(problemsubtype,"*",err,error))// &
            & " is not valid for a Helmholtz equation type of a classical field problem class."
          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_classical_field_class,problem_helmholtz_equation_type,problemsubtype]
      ELSE
        CALL flagerror("Helmholtz equation problem specification must have three entries.",err,error,*999)
      END IF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    END IF

    CALL exits("Helmholtz_ProblemSpecificationSet")
    RETURN
999 CALL errors("Helmholtz_ProblemSpecificationSet",err,error)
    CALL exits("Helmholtz_ProblemSpecificationSet")
    RETURN 1
    
  END SUBROUTINE helmholtz_problemspecificationset

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

  SUBROUTINE helmholtz_equation_problem_standard_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("HELMHOLTZ_EQUATION_PROBLEM_STANDARD_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 Helmholtz equation problem.",err,error,*999)
      END IF
      IF(problem%SPECIFICATION(3)==problem_standard_helmholtz_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 standard Helmholtz 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 standard Helmholtz 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)
            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 standard Helmholtz 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 standard Helmholtz 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 standard Helmholtz 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 standard Helmholtz equation subtype."
        CALL flagerror(local_error,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Problem is not associated.",err,error,*999)
    ENDIF
       
    exits("HELMHOLTZ_EQUATION_PROBLEM_STANDARD_SETUP")
    RETURN
999 errorsexits("HELMHOLTZ_EQUATION_PROBLEM_STANDARD_SETUP",err,error)
    RETURN 1
  END SUBROUTINE helmholtz_equation_problem_standard_setup

  !
  !================================================================================================================================
  !
 
END MODULE helmholtz_equations_routines