interface_routines.f90

Go to the documentation of this file.

MODULE interface_routines

  USE base_routines
  USE data_point_routines
  USE data_projection_routines
  USE field_routines
  USE generated_mesh_routines
  USE input_output
  USE interface_conditions_routines
  USE iso_varying_string
  USE kinds
  USE lists
  USE mesh_routines
  USE node_routines
  USE strings
  USE types

#include "macros.h"  

  IMPLICIT NONE

  PRIVATE

  !Module parameters

  !Module types

  !Module variables

  !Interfaces

  INTERFACE interface_label_get
    MODULE PROCEDURE interface_label_get_c
    MODULE PROCEDURE interface_label_get_vs
  END INTERFACE !INTERFACE_LABEL_GET
  
  INTERFACE interface_label_set
    MODULE PROCEDURE interface_label_set_c
    MODULE PROCEDURE interface_label_set_vs
  END INTERFACE !INTERFACE_LABEL_SET
  
  PUBLIC interface_mesh_add

  PUBLIC interface_create_start, interface_create_finish
  
  PUBLIC interface_coordinate_system_set,interface_coordinate_system_get
  
  PUBLIC interface_data_points_get

  PUBLIC interface_destroy, interface_mesh_connectivity_destroy, interfacepointsconnectivity_destroy

  PUBLIC interface_label_get,interface_label_set

  PUBLIC interface_nodes_get

  PUBLIC interface_mesh_connectivity_create_start, interface_mesh_connectivity_create_finish

  PUBLIC interface_user_number_find

  PUBLIC interfaces_finalise,interfaces_initialise

  PUBLIC interface_mesh_connectivity_element_xi_set,interfacemeshconnectivity_elementnumberset
  
  PUBLIC interface_mesh_connectivity_node_number_set

  PUBLIC interface_mesh_connectivity_basis_set
  
  PUBLIC interfacepointsconnectivity_createstart,interfacepointsconnectivity_createfinish
  
  PUBLIC interfacepointsconnectivity_datareprojection
  
  PUBLIC interfacepointsconnectivity_elementnumberget,interfacepointsconnectivity_elementnumberset
  
  PUBLIC interfacepointsconnectivity_pointxiget,interfacepointsconnectivity_pointxiset
  
  PUBLIC interfacepointsconnectivity_updatefromprojection
  
CONTAINS

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

  SUBROUTINE interface_mesh_add(INTERFACE,MESH,MESH_INDEX,ERR,ERROR,*)   

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(mesh_type), POINTER :: mesh
    INTEGER(INTG), INTENT(OUT) :: mesh_index
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: mesh_idx
    LOGICAL :: mesh_already_coupled
    TYPE(mesh_type), POINTER :: coupled_mesh
    TYPE(mesh_ptr_type), POINTER :: new_coupled_meshes(:)
    TYPE(region_type), POINTER :: coupled_mesh_region,mesh_region
    TYPE(varying_string) :: local_error

    NULLIFY(new_coupled_meshes)
    
    enters("INTERFACE_MESH_ADD",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        CALL flagerror("Interface has already been finished.",err,error,*999)
      ELSE
        IF(ASSOCIATED(mesh)) THEN
          IF(mesh%MESH_FINISHED) THEN
            mesh_region=>mesh%REGION
            IF(ASSOCIATED(mesh_region)) THEN
              ALLOCATE(new_coupled_meshes(interface%NUMBER_OF_COUPLED_MESHES+1),stat=err)
              IF(err/=0) CALL flagerror("Could not allocate new coupled meshes.",err,error,*999)
              !Check that the mesh is not already in the list of meshes for the interface.
              IF(interface%NUMBER_OF_COUPLED_MESHES>0) THEN
                IF(ASSOCIATED(interface%COUPLED_MESHES)) THEN
                  mesh_already_coupled=.false.
                  DO mesh_idx=1,interface%NUMBER_OF_COUPLED_MESHES
                    coupled_mesh=>interface%COUPLED_MESHES(mesh_idx)%PTR
                    IF(ASSOCIATED(coupled_mesh)) THEN
                      coupled_mesh_region=>coupled_mesh%REGION
                      IF(ASSOCIATED(coupled_mesh_region)) THEN
                        IF(mesh_region%USER_NUMBER==coupled_mesh_region%USER_NUMBER) THEN
                          IF(mesh%USER_NUMBER==coupled_mesh%USER_NUMBER) THEN
                            mesh_already_coupled=.true.
                            EXIT
                          ENDIF
                        ENDIF
                      ELSE
                        local_error="Coupled interface mesh region for mesh index "// &
                          & trim(number_to_vstring(mesh_idx,"*",err,error))//" is not associated."
                        CALL flagerror(local_error,err,error,*999)
                      ENDIF
                    ELSE
                      local_error="Coupled interface mesh for mesh index "//trim(number_to_vstring(mesh_idx,"*",err,error))// &
                        & " is not associated."
                      CALL flagerror(local_error,err,error,*999)
                    ENDIF
                    new_coupled_meshes(mesh_idx)%PTR=>interface%COUPLED_MESHES(mesh_idx)%PTR
                  ENDDO !mesh_idx
                  IF(mesh_already_coupled) THEN
                    local_error="The supplied mesh has already been added to the list of coupled meshes at mesh index "// &
                      & trim(number_to_vstring(mesh_idx,"*",err,error))//"."
                    CALL flagerror(local_error,err,error,*999)
                  ENDIF
                  DEALLOCATE(interface%COUPLED_MESHES)
                ELSE
                  CALL flagerror("Interface coupled meshes is not associated.",err,error,*999)
                ENDIF
              ENDIF
              !Add the mesh to the list of coupled meshes
              new_coupled_meshes(interface%NUMBER_OF_COUPLED_MESHES+1)%PTR=>mesh
              interface%COUPLED_MESHES=>new_coupled_meshes
              !Increment the number of coupled meshes and return the index
              interface%NUMBER_OF_COUPLED_MESHES=interface%NUMBER_OF_COUPLED_MESHES+1
              mesh_index=interface%NUMBER_OF_COUPLED_MESHES
            ELSE
              CALL flagerror("Mesh region is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Mesh has not been finished.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Mesh is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_MESH_ADD")
    RETURN
999 IF(ASSOCIATED(new_coupled_meshes)) DEALLOCATE(new_coupled_meshes)
    errorsexits("INTERFACE_MESH_ADD",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_add

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

  SUBROUTINE interface_create_finish(INTERFACE,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: local_error
     
    enters("INTERFACE_CREATE_FINISH",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        CALL flagerror("Interface has already been finished.",err,error,*999)
      ELSE
        IF(interface%NUMBER_OF_COUPLED_MESHES<2) THEN
          local_error="Invalid mesh coupling. Only "//trim(number_to_vstring(interface%NUMBER_OF_COUPLED_MESHES,"*",err,error))// &
            & " have been coupled. The number of coupled meshes must be >= 2."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
        interface%INTERFACE_FINISHED=.true.
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    IF(diagnostics1) THEN
      CALL write_string(diagnostic_output_type,"Interface :",err,error,*999)
      CALL write_string_value(diagnostic_output_type,"  User number = ",interface%USER_NUMBER,err,error,*999)
      CALL write_string_value(diagnostic_output_type,"  Global number = ",interface%GLOBAL_NUMBER,err,error,*999)
      CALL write_string_value(diagnostic_output_type,"  Label = ",interface%LABEL,err,error,*999)
      IF(ASSOCIATED(interface%INTERFACES)) THEN
        IF(ASSOCIATED(interface%INTERFACES%PARENT_REGION)) THEN
          CALL write_string_value(diagnostic_output_type,"  Parent region user number = ",interface%INTERFACES% &
            & parent_region%USER_NUMBER,err,error,*999)
          CALL write_string_value(diagnostic_output_type,"  Parent region label = ",interface%INTERFACES% &
            & parent_region%LABEL,err,error,*999)        
        ELSE
          CALL flagerror("Interfaces parent region is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface interfaces is not associated.",err,error,*999)
      ENDIF
    ENDIF
    
    exits("INTERFACE_CREATE_FINISH")
    RETURN
999 errorsexits("INTERFACE_CREATE_FINISH",err,error)
    RETURN 1
  END SUBROUTINE interface_create_finish

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

  SUBROUTINE interface_create_start(USER_NUMBER,PARENT_REGION,INTERFACE,ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: user_number
    TYPE(region_type), POINTER :: parent_region
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interface_idx
    TYPE(interface_type), POINTER :: new_interface
    TYPE(interface_ptr_type), POINTER :: new_interfaces(:)
    TYPE(varying_string) :: local_error,local_string

    NULLIFY(new_interface)
    NULLIFY(new_interfaces)
    
    enters("INTERFACE_CREATE_START",err,error,*998)

    IF(ASSOCIATED(parent_region)) THEN
      IF(ASSOCIATED(interface)) THEN
        CALL flagerror("Interface is already associated.",err,error,*998)
      ELSE
        NULLIFY(interface)
        CALL interface_user_number_find(user_number,parent_region,interface,err,error,*998)
        IF(ASSOCIATED(interface)) THEN
          local_error="Interface number "//trim(number_to_vstring(user_number,"*",err,error))// &
            & " has already been created on region number "//trim(number_to_vstring(parent_region%USER_NUMBER,"*",err,error))//"."
          CALL flagerror(local_error,err,error,*998)
        ELSE        
          NULLIFY(interface)
          !Allocate and set default interface properties.
          CALL interface_initialise(new_interface,err,error,*999)
          new_interface%USER_NUMBER=user_number
          new_interface%GLOBAL_NUMBER=parent_region%INTERFACES%NUMBER_OF_INTERFACES+1
          local_string="Interface_"//number_to_vstring(user_number,"*",err,error)
          new_interface%LABEL=char(local_string)
          IF(err/=0) GOTO 999
          new_interface%INTERFACES=>parent_region%INTERFACES
          new_interface%PARENT_REGION=>parent_region
          !Add new initerface into list of interfaces in the parent region
          ALLOCATE(new_interfaces(parent_region%INTERFACES%NUMBER_OF_INTERFACES+1),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate new interfaces.",err,error,*999)
          DO interface_idx=1,parent_region%INTERFACES%NUMBER_OF_INTERFACES
            new_interfaces(interface_idx)%PTR=>parent_region%INTERFACES%INTERFACES(interface_idx)%PTR
          ENDDO !interface_idx
          new_interfaces(parent_region%INTERFACES%NUMBER_OF_INTERFACES+1)%PTR=>new_interface
          IF(ASSOCIATED(parent_region%INTERFACES%INTERFACES)) DEALLOCATE(parent_region%INTERFACES%INTERFACES)
          parent_region%INTERFACES%INTERFACES=>new_interfaces
          parent_region%INTERFACES%NUMBER_OF_INTERFACES=parent_region%INTERFACES%NUMBER_OF_INTERFACES+1
          interface=>new_interface
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Parent region is not associated.",err,error,*998)
    ENDIF
    
    exits("INTERFACE_CREATE_START")
    RETURN
999 IF(ASSOCIATED(new_interfaces)) DEALLOCATE(new_interfaces)
    CALL interface_finalise(interface,err,error,*998)
998 errorsexits("INTERFACE_CREATE_START",err,error)
    RETURN 1
  END SUBROUTINE interface_create_start

  !
  !================================================================================================================================
  !
  
  SUBROUTINE interface_coordinate_system_get(INTERFACE,COORDINATE_SYSTEM,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(coordinate_system_type), POINTER :: coordinate_system
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    
    enters("INTERFACE_COORDINATE_SYSTEM_GET",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        IF(ASSOCIATED(coordinate_system)) THEN
          CALL flagerror("Coordinate system is already associated.",err,error,*999)
        ELSE
          coordinate_system=>interface%COORDINATE_SYSTEM
        ENDIF
      ELSE
        CALL flagerror("Interface has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_COORDINATE_SYSTEM_GET")
    RETURN
999 errorsexits("INTERFACE_COORDINATE_SYSTEM_GET",err,error)
    RETURN 1
  END SUBROUTINE interface_coordinate_system_get
  
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interface_coordinate_system_set(INTERFACE,COORDINATE_SYSTEM,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(coordinate_system_type), POINTER :: coordinate_system
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        CALL flagerror("Interface has been finished.",err,error,*999)
      ELSE
        IF(ASSOCIATED(coordinate_system)) THEN
          IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
            interface%COORDINATE_SYSTEM=>coordinate_system
          ELSE
            CALL flagerror("Coordinate system has not been finished.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Coordinate system is not associated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_COORDINATE_SYSTEM_SET")
    RETURN
999 errorsexits("INTERFACE_COORDINATE_SYSTEM_SET",err,error)
    RETURN 1
  END SUBROUTINE interface_coordinate_system_set
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interface_data_points_get(INTERFACE,DATA_POINTS,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(data_points_type), POINTER :: data_points
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
 
    enters("REGION_DATA_POINTS_GET",err,error,*998)

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN 
        IF(ASSOCIATED(data_points)) THEN
          CALL flagerror("Data points is already associated.",err,error,*998)
        ELSE
          data_points=>interface%DATA_POINTS
          IF(.NOT.ASSOCIATED(data_points)) CALL flagerror("Data points is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface has not been finished.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*998)
    ENDIF
       
    exits("INTERFACE_DATA_POINTS_GET")
    RETURN
999 NULLIFY(data_points)
998 errorsexits("INTERFACE_DATA_POINTS_GET",err,error)
    RETURN 1
    
  END SUBROUTINE interface_data_points_get    

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

  SUBROUTINE interface_destroy(INTERFACE,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interface_idx,interface_position
    TYPE(interface_ptr_type), POINTER :: new_interfaces(:)
    TYPE(interfaces_type), POINTER :: interfaces
     
    NULLIFY(new_interfaces)

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

    IF(ASSOCIATED(interface)) THEN
      interfaces=>interface%INTERFACES
      IF(ASSOCIATED(interfaces)) THEN
        interface_position=interface%GLOBAL_NUMBER

        !Destroy all the interface condition components
        CALL interface_finalise(interface,err,error,*999)
        
        !Remove the interface from the list of interfaces
        IF(interfaces%NUMBER_OF_INTERFACES>1) THEN
          ALLOCATE(new_interfaces(interfaces%NUMBER_OF_INTERFACES-1),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate new interface conditions.",err,error,*999)
          DO interface_idx=1,interfaces%NUMBER_OF_INTERFACES
            IF(interface_idx<interface_position) THEN
              new_interfaces(interface_idx)%PTR=>interfaces%INTERFACES(interface_idx)%PTR
            ELSE IF(interface_idx>interface_position) THEN
              interfaces%INTERFACES(interface_idx)%PTR%GLOBAL_NUMBER=interfaces%INTERFACES(interface_idx)%PTR%GLOBAL_NUMBER-1
              new_interfaces(interface_idx-1)%PTR=>interfaces%INTERFACES(interface_idx)%PTR
            ENDIF
          ENDDO !interface_idx
          IF(ASSOCIATED(interfaces%INTERFACES)) DEALLOCATE(interfaces%INTERFACES)
          interfaces%INTERFACES=>new_interfaces
          interfaces%NUMBER_OF_INTERFACES=interfaces%NUMBER_OF_INTERFACES-1
        ELSE
          DEALLOCATE(interfaces%INTERFACES)
          interfaces%NUMBER_OF_INTERFACES=0
        ENDIF
        
      ELSE
        CALL flagerror("Interface interfaces is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF    
    
    exits("INTERFACE_DESTROY")
    RETURN
999 errorsexits("INTERFACE_DESTROY",err,error)
    RETURN 1
  END SUBROUTINE interface_destroy

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

  SUBROUTINE interface_finalise(INTERFACE,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
     
    enters("INTERFACE_FINALISE",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(ASSOCIATED(interface%COUPLED_MESHES)) DEALLOCATE(interface%COUPLED_MESHES)
      IF(ASSOCIATED(interface%MESH_CONNECTIVITY)) &
         & CALL interface_mesh_connectivity_finalise(interface%MESH_CONNECTIVITY,err,error,*999)
      IF(ASSOCIATED(interface%pointsConnectivity)) &
        & CALL interfacepointsconnectivity_finalise(interface%pointsConnectivity,err,error,*999)
      IF(ASSOCIATED(interface%NODES)) CALL nodes_destroy(interface%NODES,err,error,*999)
      CALL meshes_finalise(interface%MESHES,err,error,*999)
      CALL fields_finalise(interface%FIELDS,err,error,*999)
      CALL interface_conditions_finalise(interface%INTERFACE_CONDITIONS,err,error,*999)
      DEALLOCATE(interface)
    ENDIF
    
    exits("INTERFACE_FINALISE")
    RETURN
999 errorsexits("INTERFACE_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE interface_finalise

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

  SUBROUTINE interface_initialise(INTERFACE,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
     
    enters("INTERFACE_INITIALISE",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      CALL flagerror("Interface is already associated.",err,error,*999)
    ELSE
      ALLOCATE(interface,stat=err)
      IF(err/=0) CALL flagerror("Could not allocate interface.",err,error,*999)
      interface%USER_NUMBER=0
      interface%GLOBAL_NUMBER=0
      interface%INTERFACE_FINISHED=.false.
      interface%LABEL=""
      NULLIFY(interface%INTERFACES)
      NULLIFY(interface%PARENT_REGION)
      interface%NUMBER_OF_COUPLED_MESHES=0
      NULLIFY(interface%COUPLED_MESHES)
      NULLIFY(interface%MESH_CONNECTIVITY)
      NULLIFY(interface%pointsConnectivity)
      NULLIFY(interface%NODES)
      NULLIFY(interface%MESHES)
      NULLIFY(interface%GENERATED_MESHES)
      NULLIFY(interface%FIELDS)
      NULLIFY(interface%INTERFACE_CONDITIONS)
      NULLIFY(interface%COORDINATE_SYSTEM)
      NULLIFY(interface%DATA_POINTS)
      CALL meshes_initialise(interface,err,error,*999)
      CALL generated_meshes_initialise(interface,err,error,*999)
      CALL fields_initialise(interface,err,error,*999)
      CALL interface_conditions_initialise(interface,err,error,*999)
    ENDIF
    
    exits("INTERFACE_INITIALISE")
    RETURN
999 errorsexits("INTERFACE_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE interface_initialise

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

  SUBROUTINE interface_label_get_c(INTERFACE,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    CHARACTER(LEN=*), INTENT(OUT) :: label
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: c_length,vs_length

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

    IF(ASSOCIATED(interface)) THEN
      c_length=len(label)
      vs_length=len_trim(interface%LABEL)
      IF(c_length>vs_length) THEN
        label=char(interface%LABEL,vs_length)
      ELSE
        label=char(interface%LABEL,c_length)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_LABEL_GET_C")
    RETURN
999 errorsexits("INTERFACE_LABEL_GET_C",err,error)
    RETURN 1
    
  END SUBROUTINE interface_label_get_c

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

  SUBROUTINE interface_label_get_vs(INTERFACE,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(varying_string), INTENT(OUT) :: label
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      !\todo The following line crashes the AIX compiler unless it has a VAR_STR(CHAR()) around it
      label=var_str(char(interface%LABEL))
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_LABEL_GET_VS")
    RETURN
999 errorsexits("INTERFACE_LABEL_GET_VS",err,error)
    RETURN 1
    
  END SUBROUTINE interface_label_get_vs

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

  SUBROUTINE interface_label_set_c(INTERFACE,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    CHARACTER(LEN=*), INTENT(IN) :: label
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        CALL flagerror("Interface has been finished.",err,error,*999)
      ELSE
        interface%LABEL=label
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_LABEL_SET_C")
    RETURN
999 errorsexits("INTERFACE_LABEL_SET_C",err,error)
    RETURN 1
  END SUBROUTINE interface_label_set_c

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

  SUBROUTINE interface_label_set_vs(INTERFACE,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(varying_string), INTENT(IN) :: label
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        CALL flagerror("Interface has been finished.",err,error,*999)
      ELSE
        interface%LABEL=label
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_LABEL_SET_VS")
    RETURN
999 errorsexits("INTERFACE_LABEL_SET_VS",err,error)
    RETURN 1
  END SUBROUTINE interface_label_set_vs

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

  SUBROUTINE interface_nodes_get(INTERFACE,NODES,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(nodes_type), POINTER :: nodes
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
 
    enters("INTERFACE_NODES_GET",err,error,*998)

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN 
        IF(ASSOCIATED(nodes)) THEN
          CALL flagerror("Nodes is already associated.",err,error,*998)
        ELSE
          nodes=>interface%NODES
          IF(.NOT.ASSOCIATED(nodes)) CALL flagerror("Nodes is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface has not been finished.",err,error,*998)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*998)
    ENDIF
       
    exits("INTERFACE_NODES_GET")
    RETURN
999 NULLIFY(nodes)
998 errorsexits("INTERFACE_NODES_GET",err,error)
    RETURN 1
    
  END SUBROUTINE interface_nodes_get

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

  SUBROUTINE interface_mesh_connectivity_create_finish(INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interfaceelementidx,coupledmeshidx
    TYPE(varying_string) :: local_error

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

    IF(ASSOCIATED(interface_mesh_connectivity)) THEN
      IF(interface_mesh_connectivity%MESH_CONNECTIVITY_FINISHED) THEN
        CALL flagerror("Interface meshes connectivity has already been finished.",err,error,*999)
      ELSE
        !Check if connectivity from each interface element to an appropriate element in the coupled meshes has been setup
        DO interfaceelementidx=1,interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS
          DO coupledmeshidx=1,interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES
            IF (interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)% &
              & coupled_mesh_element_number==0) THEN
              local_error="The connectivity from interface element " &
                //trim(number_to_vstring(interfaceelementidx,"*",err,error))//" to an element in coupled mesh " & 
                //trim(number_to_vstring(coupledmeshidx,"*",err,error))//" has not been defined."
              CALL flagerror(local_error,err,error,*999)
            ENDIF
          ENDDO !CoupledMeshIdx
        ENDDO !InterfaceElementIdx
        !Calculate line or face numbers for coupled mesh elements that are connected to the interface mesh
        CALL interfacemeshconnectivity_connectedlinescalculate(interface_mesh_connectivity,err,error,*999)
        interface_mesh_connectivity%MESH_CONNECTIVITY_FINISHED=.true.
      ENDIF
    ELSE
      CALL flagerror("Interface meshes connectivity is not associated.",err,error,*999)
    ENDIF

    exits("INTERFACE_MESH_CONNECTIVITY_CREATE_FINISH")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_CREATE_FINISH",err,error)
    RETURN 1
  END SUBROUTINE interface_mesh_connectivity_create_finish

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

  SUBROUTINE interface_mesh_connectivity_create_start(INTERFACE,MESH,INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(mesh_type), POINTER :: mesh
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        IF(ASSOCIATED(interface%MESH_CONNECTIVITY)) THEN
          CALL flagerror("The interface already has a meshes connectivity associated.",err,error,*999)
        ELSE
          !Initialise the meshes connectivity
          CALL interface_mesh_connectivity_initialise(interface,mesh,err,error,*999)
          !Return the pointer
          interface_mesh_connectivity=>interface%MESH_CONNECTIVITY
        ENDIF
      ELSE
        CALL flagerror("Interface has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_MESH_CONNECTIVITY_CREATE_START")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_CREATE_START",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_connectivity_create_start

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

  SUBROUTINE interface_mesh_connectivity_destroy(INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    
    enters("INTERFACE_MESH_CONNECTIVITY_DESTROY",err,error,*999)

    IF(ASSOCIATED(interface_mesh_connectivity)) THEN
      CALL interface_mesh_connectivity_finalise(interface_mesh_connectivity,err,error,*999)
    ELSE
      CALL flagerror("Interface meshes connectivity is not associated.",err,error,*999)
    ENDIF
       
    exits("INTERFACE_MESH_CONNECTIVITY_DESTROY")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_DESTROY",err,error)
    RETURN 1
  END SUBROUTINE interface_mesh_connectivity_destroy

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

  SUBROUTINE interface_mesh_connectivity_basis_set(INTERFACE_MESH_CONNECTIVITY,BASIS,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    TYPE(basis_type), POINTER :: basis
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(interface_element_connectivity_type), POINTER :: element_connectivity
    INTEGER(INTG) :: interfaceelementidx,coupledmeshidx,numberofinterfaceelementnodes,numberofcoupledmeshxidirections

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

    IF(ASSOCIATED(interface_mesh_connectivity)) THEN
      IF(interface_mesh_connectivity%MESH_CONNECTIVITY_FINISHED) THEN
        CALL flagerror("Interface mesh connectivity already been finished.",err,error,*999)
      ELSE
        IF(ASSOCIATED(interface_mesh_connectivity%BASIS)) THEN
          CALL flagerror("Mesh connectivity basis already associated.",err,error,*999)
        ELSE
          IF(ASSOCIATED(basis)) THEN
            interface_mesh_connectivity%BASIS=>basis
            !Now that the mesh connectivity basis is set the number of interface element nodes can be determined and now MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(InterfaceElementIdx,CoupledMeshIdx)%XI can be allocated
            !\todo NumberOfCoupledMeshXiDirections currently set to the number of interface mesh xi directions + 1. Restructure ELEMENT_CONNECTIVITY type see below
            numberofcoupledmeshxidirections=interface_mesh_connectivity%INTERFACE_MESH%NUMBER_OF_DIMENSIONS+1
            numberofinterfaceelementnodes=interface_mesh_connectivity%BASIS%NUMBER_OF_NODES
            DO interfaceelementidx=1,interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS
              DO coupledmeshidx=1,interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES
                element_connectivity=>interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)
                !\todo Update mesh component index to look at the number of mesh components in each element. 
                !\todo Currently this defaults to the first mesh component ie %XI(NumberOfInterfaceMeshXi,1,NumberOfInterfaceElementNodes)). 
                !\todo The interface mesh types will also need to be restructured.
                !eg. INTERFACE%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(InterfaceElementIdx,CoupledMeshIdx)%MESH_COMPONENT(MeshComponentIdx)%XI(NumberOfCoupledMeshXiDirections,NumberOfInterfaceElementNodes) and adding appropriate initialize/finialize routines
                ALLOCATE(element_connectivity%XI(numberofcoupledmeshxidirections,1,numberofinterfaceelementnodes),stat=err)
                IF(err/=0) CALL flagerror("Could not allocate interface element connectivity.",err,error,*999)
                element_connectivity%XI=0.0_dp
              ENDDO !CoupledMeshIdx
            ENDDO !InterfaceElementIdx
          ELSE
            CALL flagerror("Basis to set mesh connectivity not associated.",err,error,*999)
          ENDIF
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Interface mesh connectivity is not associated.",err,error,*999)
    ENDIF

    exits("INTERFACE_MESH_CONNECTIVITY_BASIS_SET")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_BASIS_SET",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_connectivity_basis_set

  !
  !================================================================================================================================
  !
    
  SUBROUTINE interface_mesh_connectivity_element_xi_set(INTERFACE_MESH_CONNECTIVITY,INTERFACE_MESH_ELEMENT_NUMBER, &
    & coupled_mesh_index,coupled_mesh_element_number,interface_mesh_local_node_number,interface_mesh_component_number,xi, &
    & err,error,*)

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(IN) :: interface_mesh_element_number
    INTEGER(INTG), INTENT(IN) :: coupled_mesh_index
    INTEGER(INTG), INTENT(IN) :: coupled_mesh_element_number
    INTEGER(INTG), INTENT(IN) :: interface_mesh_local_node_number
    INTEGER(INTG), INTENT(IN) :: interface_mesh_component_number
    REAL(DP), INTENT(IN) :: xi(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: coupled_mesh_number_of_xi
    TYPE(interface_element_connectivity_type), POINTER :: element_connectivity
    
    enters("INTERFACE_MESH_CONNECTIVITY_ELEMENT_XI_SET",err,error,*999)

    IF(ASSOCIATED(interface_mesh_connectivity)) THEN
      IF(interface_mesh_connectivity%MESH_CONNECTIVITY_FINISHED) THEN
        CALL flagerror("Interface mesh connectivity already been finished.",err,error,*999)
      ELSE
        IF(ALLOCATED(interface_mesh_connectivity%ELEMENT_CONNECTIVITY)) THEN
          IF((interface_mesh_element_number>0).AND. &
            & (interface_mesh_element_number<=interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS)) THEN
            IF((coupled_mesh_index>0).AND.(coupled_mesh_index<=interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES)) THEN
              IF((coupled_mesh_element_number>0).AND.(coupled_mesh_element_number<= &
                & interface_mesh_connectivity%INTERFACE%COUPLED_MESHES(coupled_mesh_index)%PTR%NUMBER_OF_ELEMENTS))THEN
                IF((interface_mesh_component_number>0).AND. &
                  & (interface_mesh_component_number<=interface_mesh_connectivity%INTERFACE_MESH%NUMBER_OF_COMPONENTS)) THEN
                  IF((interface_mesh_local_node_number>0).AND.(interface_mesh_local_node_number<= &
                    & interface_mesh_connectivity%BASIS%NUMBER_OF_NODES))THEN
                    element_connectivity=>interface_mesh_connectivity% &
                      & element_connectivity(interface_mesh_element_number,coupled_mesh_index)
                    IF(element_connectivity%COUPLED_MESH_ELEMENT_NUMBER==coupled_mesh_element_number)THEN
                      coupled_mesh_number_of_xi = interface_mesh_connectivity%INTERFACE%COUPLED_MESHES(coupled_mesh_index)%PTR% & 
                        & topology(1)%PTR%ELEMENTS%ELEMENTS(coupled_mesh_element_number)%BASIS%NUMBER_OF_XI
                      IF(SIZE(xi)==coupled_mesh_number_of_xi) THEN
                        !\todo the ELEMENT_CONNECTIVITY%XI array needs to be restructured to efficiently utilize memory when coupling bodies with 2xi directions to bodies with 3xi directions using an interface.
                        element_connectivity%XI(1:coupled_mesh_number_of_xi,interface_mesh_component_number, &
                          & interface_mesh_local_node_number)=xi(1:coupled_mesh_number_of_xi)
                      ELSE
                        CALL flagerror("The size of the xi array provided does not match the coupled mesh element's' number"// &
                          & " of xi.",err,error,*999)
                      ENDIF
                    ELSE
                      CALL flagerror("Coupled mesh element number doesn't match that set to the interface.",err,error,*999)
                    ENDIF
                  ELSE
                    CALL flagerror("Interface local node number is out of range.",err,error,*999)
                  ENDIF
                ELSE
                  CALL flagerror("Interface component number is out of range.",err,error,*999)
                ENDIF
              ELSE
                CALL flagerror("Coupled mesh element number out of range.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Interface coupled mesh index number out of range.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Interface mesh element number out of range.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Interface elements connectivity array not allocated.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Interface mesh connectivity is not associated.",err,error,*999)
    ENDIF

    exits("INTERFACE_MESH_CONNECTIVITY_ELEMENT_XI_SET")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_ELEMENT_XI_SET",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_connectivity_element_xi_set

  !
  !================================================================================================================================
  !
  
  SUBROUTINE interfacemeshconnectivity_elementnumberset(INTERFACE_MESH_CONNECTIVITY,INTERFACE_MESH_ELEMENT_NUMBER, &
      & coupled_mesh_index,coupled_mesh_element_number,err,error,*)

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(IN) :: interface_mesh_element_number
    INTEGER(INTG), INTENT(IN) :: coupled_mesh_index
    INTEGER(INTG), INTENT(IN) :: coupled_mesh_element_number
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(interface_element_connectivity_type), POINTER :: element_connectivity
    
    enters("InterfaceMeshConnectivity_ElementNumberSet",err,error,*999)

    IF(ASSOCIATED(interface_mesh_connectivity)) THEN
      IF(interface_mesh_connectivity%MESH_CONNECTIVITY_FINISHED) THEN
        CALL flagerror("Interface mesh connectivity has already been finished.",err,error,*999)
      ELSE
        IF((interface_mesh_element_number>0).AND.(interface_mesh_element_number<= &
          & interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS)) THEN
          IF((coupled_mesh_index>0).AND.(coupled_mesh_index<=interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES)) THEN
            IF (ALLOCATED(interface_mesh_connectivity%ELEMENT_CONNECTIVITY)) THEN
              element_connectivity=>interface_mesh_connectivity% &
                & element_connectivity(interface_mesh_element_number,coupled_mesh_index)
              element_connectivity%COUPLED_MESH_ELEMENT_NUMBER=coupled_mesh_element_number
            ELSE
              CALL flagerror("Interface elements connectivity array not allocated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Interface coupled mesh index number out of range.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Interface mesh element number out of range.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Interface mesh connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfaceMeshConnectivity_ElementNumberSet")
    RETURN
999 errorsexits("InterfaceMeshConnectivity_ElementNumberSet",err,error)
    RETURN 1
    
  END SUBROUTINE interfacemeshconnectivity_elementnumberset

  !
  !================================================================================================================================
  !
  
  SUBROUTINE interface_mesh_connectivity_node_number_set(NODES,INTERFACE_MESH_NODE_NUMBERS, &
      & first_coupled_mesh_index,first_coupled_mesh_node_numbers,second_coupled_mesh_index,second_coupled_mesh_node_numbers, &
      & err,error,*)        
        
    !Argument variables
    TYPE(nodes_type), POINTER :: nodes
    INTEGER(INTG), INTENT(IN) :: interface_mesh_node_numbers(:)
    INTEGER(INTG), INTENT(IN) :: first_coupled_mesh_index,second_coupled_mesh_index
    INTEGER(INTG), INTENT(IN) :: first_coupled_mesh_node_numbers(:),second_coupled_mesh_node_numbers(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: nodeindex
    
    enters("INTERFACE_MESH_CONNECTIVITY_NODE_NUMBER_SET",err,error,*999)

    IF(ASSOCIATED(nodes)) THEN
      IF(nodes%INTERFACE%MESH_CONNECTIVITY%MESH_CONNECTIVITY_FINISHED) THEN
        print *, 'CHECK how to circumvent! interface_routines.f90:1133'
        CALL flagerror("Interface mesh connectivity has already been finished.",err,error,*999)
      ELSE
        !Default to two coupled meshes
        ALLOCATE(nodes%COUPLED_NODES(2,SIZE(interface_mesh_node_numbers(:))))
        DO nodeindex=1,SIZE(interface_mesh_node_numbers(:))
          nodes%COUPLED_NODES(first_coupled_mesh_index,nodeindex)=first_coupled_mesh_node_numbers(nodeindex)
          nodes%COUPLED_NODES(second_coupled_mesh_index,nodeindex)=second_coupled_mesh_node_numbers(nodeindex)
        ENDDO
      ENDIF
    ELSE
      CALL flagerror("Nodes are not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_MESH_CONNECTIVITY_NODE_NUMBER_SET")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_NODE_NUMBER_SET",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_connectivity_node_number_set

  !
  !================================================================================================================================
  !
  
  SUBROUTINE interfacemeshconnectivity_connectedlinescalculate(INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_mesh_connectivity_type), POINTER :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local variables
    REAL(DP) :: xidifference
    INTEGER(INTG) :: coupledmeshidx,interfaceelementidx,coupledmeshxiidx,numberofinterfaceelementnodes,numberofinterfacemeshxi
    TYPE(interface_element_connectivity_type), POINTER :: element_connectivity
    TYPE(varying_string) :: local_error
    
    enters("InterfaceMeshConnectivity_ConnectedLinesCalculate",err,error,*999)
    
    DO coupledmeshidx=1,interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES
      DO interfaceelementidx=1,interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS
        element_connectivity=>interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)
        numberofinterfaceelementnodes=interface_mesh_connectivity%BASIS%NUMBER_OF_NODES
        numberofinterfacemeshxi=interface_mesh_connectivity%INTERFACE_MESH%NUMBER_OF_DIMENSIONS
        SELECTCASE(numberofinterfacemeshxi)
        CASE(1) !Lines
          DO coupledmeshxiidx=1,interface_mesh_connectivity%INTERFACE%COUPLED_MESHES(coupledmeshidx)%PTR%NUMBER_OF_DIMENSIONS
            ! Calculate difference between first node and last node of an element
            xidifference=element_connectivity%XI(coupledmeshxiidx,1,1)- &
              & element_connectivity%XI(coupledmeshxiidx,1,numberofinterfaceelementnodes)
            IF(abs(xidifference)<zero_tolerance) THEN
              IF(abs(element_connectivity%XI(coupledmeshxiidx,1,numberofinterfaceelementnodes))<zero_tolerance) THEN
                element_connectivity%CONNECTED_LINE=3-(coupledmeshxiidx-1)*2
              ELSE
                element_connectivity%CONNECTED_LINE=4-(coupledmeshxiidx-1)*2
              ENDIF
            ENDIF
          ENDDO
        CASE(2) !Faces
          DO coupledmeshxiidx=1,interface_mesh_connectivity%INTERFACE%COUPLED_MESHES(coupledmeshidx)%PTR%NUMBER_OF_DIMENSIONS
            xidifference=element_connectivity%XI(coupledmeshxiidx,1,1)- &
              & element_connectivity%XI(coupledmeshxiidx,1,numberofinterfaceelementnodes)
            IF(abs(xidifference)<zero_tolerance) THEN
              IF(abs(element_connectivity%XI(coupledmeshxiidx,1,numberofinterfaceelementnodes))<zero_tolerance) THEN
                element_connectivity%CONNECTED_FACE=(coupledmeshxiidx-1)*2+2
              ELSE
                element_connectivity%CONNECTED_FACE=(coupledmeshxiidx-1)*2+1
              ENDIF
            ENDIF
          ENDDO
        CASE DEFAULT 
          local_error="Number of interface mesh dimension of "//trim(number_to_vstring(numberofinterfacemeshxi,"*",err,error))// &
            & "is invalid"
          CALL flagerror(local_error,err,error,*999)
        ENDSELECT
      ENDDO
    ENDDO
    
    exits("InterfaceMeshConnectivity_ConnectedLinesCalculate")
    RETURN
999 errors("InterfaceMeshConnectivity_ConnectedLinesCalculate",err,error)
    exits("InterfaceMeshConnectivity_ConnectedLinesCalculate")
    RETURN 1
  
  END SUBROUTINE interfacemeshconnectivity_connectedlinescalculate

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

  SUBROUTINE interface_mesh_connectivity_finalise(INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_mesh_connectivity_type) :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
     
    enters("INTERFACE_MESH_CONNECTIVITY_FINALISE",err,error,*999)

    CALL interfacemeshconnectivity_elementfinalise(interface_mesh_connectivity,err,error,*999)
    NULLIFY(interface_mesh_connectivity%INTERFACE)
    NULLIFY(interface_mesh_connectivity%INTERFACE_MESH)
    NULLIFY(interface_mesh_connectivity%BASIS)
    interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS=0
    interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES=0
       
    exits("INTERFACE_MESH_CONNECTIVITY_FINALISE")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_FINALISE",err,error)
    RETURN 1
    
  END SUBROUTINE interface_mesh_connectivity_finalise

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

  SUBROUTINE interface_mesh_connectivity_initialise(INTERFACE,MESH,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(mesh_type), POINTER :: mesh
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(ASSOCIATED(interface%MESH_CONNECTIVITY)) THEN
        CALL flagerror("Interface mesh connectivity is already associated.",err,error,*999)
      ELSE
        ALLOCATE(interface%MESH_CONNECTIVITY,stat=err)
        IF(err/=0) CALL flagerror("Could not allocate interface mesh connectivity.",err,error,*999)
        interface%MESH_CONNECTIVITY%INTERFACE=>INTERFACE
        interface%MESH_CONNECTIVITY%MESH_CONNECTIVITY_FINISHED=.false.
        interface%MESH_CONNECTIVITY%INTERFACE_MESH=>mesh
        NULLIFY(interface%MESH_CONNECTIVITY%BASIS)
        CALL interfacemeshconnectivity_elementinitialise(interface,err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_MESH_CONNECTIVITY_INITIALISE")
    RETURN
999 errorsexits("INTERFACE_MESH_CONNECTIVITY_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE interface_mesh_connectivity_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_coupledelementscalculate(interfacePointsConnectivity,coupledMeshIdx,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(IN) :: coupledmeshidx
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(list_type), POINTER :: elementnumberslist
    INTEGER(INTG) :: elementidx,datapointidx,globaldatapointnumber,globalelementnumber,numberofelementdatapoints, &
      & numberOfCoupledElements,coupledElementIdx
    INTEGER(INTG), ALLOCATABLE :: elementnumbers(:)
  
    enters("InterfacePointsConnectivity_CoupledElementsCalculate",err,error,*999)

    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      IF(ALLOCATED(interfacepointsconnectivity%coupledElements)) THEN
        interfacepointsconnectivity%maxNumberOfCoupledElements(coupledmeshidx)=0; !Initialise the number of coupled mesh elements 
        DO elementidx=1,SIZE(interfacepointsconnectivity%coupledElements,1) !Number of interface elements
          numberofelementdatapoints=interfacepointsconnectivity%interfaceMesh%TOPOLOGY(1)%PTR%dataPoints% &
            & elementdatapoint(elementidx)%numberOfProjectedData !Get the number of data points in interface mesh element
          !Set up list
          NULLIFY(elementnumberslist)
          CALL list_create_start(elementnumberslist,err,error,*999)
          CALL list_data_type_set(elementnumberslist,list_intg_type,err,error,*999)
          CALL list_initial_size_set(elementnumberslist,numberofelementdatapoints,err,error,*999)
          CALL list_create_finish(elementnumberslist,err,error,*999)
          DO datapointidx=1,numberofelementdatapoints
            globaldatapointnumber=interfacepointsconnectivity%interfaceMesh%TOPOLOGY(1)%PTR%dataPoints% &
              & elementdatapoint(elementidx)%dataIndices(datapointidx)%globalNumber
            globalelementnumber=interfacepointsconnectivity%pointsConnectivity(globaldatapointnumber,coupledmeshidx)% &
              & coupledmeshelementnumber
            CALL list_item_add(elementnumberslist,globalelementnumber,err,error,*999)
          ENDDO !dataPointIdx
          CALL list_remove_duplicates(elementnumberslist,err,error,*999)
          CALL list_detach_and_destroy(elementnumberslist,numberofcoupledelements,elementnumbers, &
            & err,error,*999)
          IF(ALLOCATED(interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%elementNumbers)) & 
            & DEALLOCATE(interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%elementNumbers) !for updating coupledElements after projection
          ALLOCATE(interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)% &
            & elementnumbers(numberofcoupledelements),stat=err)
          IF(err/=0) CALL flagerror("Could not allocate coupled mesh element numbers.",err,error,*999)
          DO coupledelementidx=1,numberofcoupledelements
            interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%elementNumbers(coupledelementidx)= &
              & elementnumbers(coupledelementidx)
          ENDDO
          interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%numberOfCoupledElements=numberofcoupledelements
          IF(ALLOCATED(elementnumbers)) DEALLOCATE(elementnumbers)
          IF(interfacepointsconnectivity%maxNumberOfCoupledElements(coupledmeshidx)<numberofcoupledelements) THEN
            interfacepointsconnectivity%maxNumberOfCoupledElements(coupledmeshidx)=numberofcoupledelements
          ENDIF
        ENDDO !elementIdx
      ELSE
        CALL flagerror("Interface points connectivity coupled elements is not allocated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_CoupledElementsCalculate")
    RETURN
999 errors("InterfacePointsConnectivity_CoupledElementsCalculate",err,error)
    exits("InterfacePointsConnectivity_CoupledElementsCalculate")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_coupledelementscalculate
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_coupledelementsfinalise(interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: elementidx,coupledmeshidx
  
    enters("InterfacePointsConnectivity_CoupledElementsFinalise",err,error,*999)

    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      IF(ALLOCATED(interfacepointsconnectivity%coupledElements)) THEN
        DO coupledmeshidx=1,SIZE(interfacepointsconnectivity%coupledElements,2)
          DO elementidx=1,SIZE(interfacepointsconnectivity%coupledElements,1)
            interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%numberOfCoupledElements=0
            IF(ALLOCATED(interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%elementNumbers)) &
              & DEALLOCATE(interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%elementNumbers)
          ENDDO
        ENDDO
        DEALLOCATE(interfacepointsconnectivity%coupledElements)
      END IF
      IF(ALLOCATED(interfacepointsconnectivity%maxNumberOfCoupledElements)) &
        & DEALLOCATE(interfacepointsconnectivity%maxNumberOfCoupledElements)
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_CoupledElementsFinalise")
    RETURN
999 errors("InterfacePointsConnectivity_CoupledElementsFinalise",err,error)
    exits("InterfacePointsConnectivity_CoupledElementsFinalise")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_coupledelementsfinalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_coupledelementsinitialise(interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: numberofinterfaceelements,numberofcoupledmeshes,coupledmeshidx,elementidx
    INTEGER(INTG) :: dummyerr
    TYPE(varying_string)  :: dummyerror
  
    enters("InterfacePointsConnectivity_CoupledElementsInitialise",err,error,*999)

     IF(ASSOCIATED(interfacepointsconnectivity)) THEN
       IF(ALLOCATED(interfacepointsconnectivity%coupledElements)) THEN
         CALL flagerror("Interface points connectivity coupled elements is already allocated.",err,error,*999)
       ELSE
         IF(ASSOCIATED(interfacepointsconnectivity%interface)) THEN
           IF(ASSOCIATED(interfacepointsconnectivity%interfaceMesh)) THEN
             numberofinterfaceelements=interfacepointsconnectivity%interfaceMesh%NUMBER_OF_ELEMENTS
             numberofcoupledmeshes=interfacepointsconnectivity%interface%NUMBER_OF_COUPLED_MESHES
             ALLOCATE(interfacepointsconnectivity%coupledElements(numberofinterfaceelements,numberofcoupledmeshes),stat=err)
             IF(err/=0) CALL flagerror("Could not allocate points connectivity coupled element.",err,error,*999)
             DO coupledmeshidx=1,numberofcoupledmeshes
               DO elementidx=1,numberofinterfaceelements
                 interfacepointsconnectivity%coupledElements(elementidx,coupledmeshidx)%numberOfCoupledElements=0
               ENDDO !elementIdx
             ENDDO !coupledMeshIdx
           ELSE
             CALL flagerror("Interface mesh is not associated.",err,error,*999)
           ENDIF
         ELSE
            CALL flagerror("Interface is not associated.",err,error,*999)
         ENDIF
       ENDIF
     ELSE
       CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
     ENDIF
    
    exits("InterfacePointsConnectivity_CoupledElementsInitialise")
    RETURN
999 CALL interfacepointsconnectivity_coupledelementsfinalise(interfacepointsconnectivity,dummyerr,dummyerror,*998) 
998 errors("InterfacePointsConnectivity_CoupledElementsInitialise",err,error)
    exits("InterfacePointsConnectivity_CoupledElementsInitialise")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_coupledelementsinitialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_createfinish(interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG) :: coupledmeshidx
  
    enters("InterfacePointsConnectivity_CreateFinish",err,error,*999)

    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      IF(interfacepointsconnectivity%pointsConnectivityFinished) THEN
        CALL flagerror("Interface points connectivity has already been finished.",err,error,*999)
      ELSE
        CALL interfacepointsconnectivity_reducedxicalculate(interfacepointsconnectivity,err,error,*999) 
        interface=>interfacepointsconnectivity%interface
        IF(ASSOCIATED(interface)) THEN
          CALL interfacepointsconnectivity_coupledelementsinitialise(interfacepointsconnectivity,err,error,*999) 
          DO coupledmeshidx=1,interfacepointsconnectivity%interface%NUMBER_OF_COUPLED_MESHES
            CALL interfacepointsconnectivity_coupledelementscalculate(interfacepointsconnectivity,coupledmeshidx,err,error,*999) 
          ENDDO  
        ELSE
          CALL flagerror("Interface is not associated.",err,error,*999)
        ENDIF
        interfacepointsconnectivity%pointsConnectivityFinished=.true.   
      ENDIF
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_CreateFinish")
    RETURN
999 errorsexits("InterfacePointsConnectivity_CreateFinish",err,error)
    RETURN 1
  END SUBROUTINE interfacepointsconnectivity_createfinish

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

  SUBROUTINE interfacepointsconnectivity_createstart(interface,interfaceMesh,interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(mesh_type), POINTER :: interfacemesh
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interface)) THEN
      IF(interface%INTERFACE_FINISHED) THEN
        IF(ASSOCIATED(interface%pointsConnectivity)) THEN
          CALL flagerror("The interface already has a points connectivity associated.",err,error,*999)
        ELSE
          IF(ASSOCIATED(interfacemesh)) THEN
            CALL interfacepointsconnectivity_initialise(interface,interfacemesh,err,error,*999)
            !Return the pointer
            interfacepointsconnectivity=>interface%pointsConnectivity
          ELSE
            CALL flagerror("Interface mesh is not associated.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Interface has not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_CreateStart")
    RETURN
999 errorsexits("InterfacePointsConnectivity_CreateStart",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_createstart
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_datareprojection(interface,interfaceCondition,err,error,*) 

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(interface_condition_type), POINTER :: interfacecondition
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    TYPE(field_type), POINTER :: dependentfieldfixed,dependentfieldprojection
    TYPE(data_points_type), POINTER :: datapoints
    TYPE(data_projection_type), POINTER :: dataprojection
    TYPE(field_interpolated_point_ptr_type), POINTER :: interpolatedpoints(:)
    TYPE(field_interpolated_point_type), POINTER :: interpolatedpoint
    TYPE(field_interpolation_parameters_ptr_type), POINTER :: interpolationparameters(:)
    INTEGER(INTG) :: fixedbodyidx,projectionbodyidx,datapointidx
    INTEGER(INTG) :: elementnumber,numberofgeometriccomponents
    INTEGER(INTG) :: coupledmeshfacelinenumber,component
  
    enters("InterfacePointsConnectivity_DataReprojection",err,error,*999)
    
    NULLIFY(interpolatedpoints)
    NULLIFY(interpolationparameters)
    fixedbodyidx=2 !\todo: need to generalise
    projectionbodyidx=1

    IF(ASSOCIATED(interface)) THEN
      IF(ASSOCIATED(interfacecondition)) THEN
        interfacepointsconnectivity=>interface%pointsConnectivity
        IF(ASSOCIATED(interfacepointsconnectivity)) THEN
          datapoints=>interface%DATA_POINTS
          IF(ASSOCIATED(datapoints)) THEN

            !Evaluate data points positions
            dependentfieldfixed=>interfacecondition%DEPENDENT%FIELD_VARIABLES(fixedbodyidx)%PTR%FIELD
            IF(ASSOCIATED(dependentfieldfixed)) THEN
              numberofgeometriccomponents=dependentfieldfixed%GEOMETRIC_FIELD%VARIABLES(1)%NUMBER_OF_COMPONENTS
              CALL field_interpolation_parameters_initialise(dependentfieldfixed,interpolationparameters,err,error,*999, &
                & field_geometric_components_type)
              CALL field_interpolated_points_initialise(interpolationparameters,interpolatedpoints,err,error,*999, &
                & field_geometric_components_type)
              interpolatedpoint=>interpolatedpoints(field_u_variable_type)%PTR
              DO datapointidx=1,datapoints%NUMBER_OF_DATA_POINTS
                elementnumber=interfacepointsconnectivity%pointsConnectivity(datapointidx,fixedbodyidx)% &
                  & coupledmeshelementnumber
                coupledmeshfacelinenumber=dependentfieldfixed%DECOMPOSITION%TOPOLOGY%ELEMENTS% &
                  & elements(elementnumber)% &
                  & element_faces(interfacepointsconnectivity%pointsConnectivity(datapointidx,fixedbodyidx)% &
                  & elementlinefacenumber)
                CALL field_interpolation_parameters_face_get(field_values_set_type,coupledmeshfacelinenumber, &
                  & interpolationparameters(field_u_variable_type)%PTR,err,error,*999,field_geometric_components_type)
                CALL field_interpolate_xi(no_part_deriv,interfacepointsconnectivity%pointsConnectivity(datapointidx, &
                  & fixedbodyidx)%reducedXi(:),interpolatedpoint,err,error,*999,field_geometric_components_type) !Interpolate contact data points on each surface
                DO component=1,numberofgeometriccomponents
                  datapoints%DATA_POINTS(datapointidx)%position(component) = interpolatedpoint%VALUES(component,no_part_deriv)
                ENDDO !component
              ENDDO !dataPointIdx
            ELSE
              CALL flagerror("Fixed dependent field is not associated.",err,error,*999)
            ENDIF
            
            !Data reprojection and update points connectivity information with the projection results
            dataprojection=>datapoints%DATA_PROJECTIONS(projectionbodyidx+1)%PTR 
            CALL write_string(general_output_type,"ProjectedBodyDataProjectionLabel",err,error,*999)
            CALL write_string(general_output_type,dataprojection%label,err,error,*999)
            IF(ASSOCIATED(dataprojection)) THEN
              dependentfieldprojection=>interfacecondition%DEPENDENT%FIELD_VARIABLES(projectionbodyidx)%PTR%FIELD
              IF(ASSOCIATED(dependentfieldprojection)) THEN
                !Projection the data points (with know spatial positions) on the projection dependent field 
                CALL dataprojection_datapointsprojectionevaluate(dataprojection,dependentfieldprojection,err,error,*999)
                CALL interfacepointsconnectivity_updatefromprojection(interfacepointsconnectivity,dataprojection, &
                  & projectionbodyidx,err,error,*999) 
              ELSE
                CALL flagerror("Projection dependent field is not associated.",err,error,*999)
              ENDIF       
            ELSE
              CALL flagerror("Interface data projection is not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Interface data points is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface condition is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_DataReprojection")
    RETURN
999 errors("InterfacePointsConnectivity_DataReprojection",err,error)
    exits("InterfacePointsConnectivity_DataReprojection")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_datareprojection
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_destroy(interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

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

    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      CALL interfacepointsconnectivity_finalise(interfacepointsconnectivity,err,error,*999) 
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_Destroy")
    RETURN
999 errorsexits("InterfacePointsConnectivity_Destroy",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_destroy
  
  !
  !================================================================================================================================
  !
  
  SUBROUTINE interfacepointsconnectivity_elementnumberget(pointsConnectivity,dataPointUserNumber,coupledMeshIndexNumber, &
      & meshcomponentnumber,coupledmeshuserelementnumber,err,error,*)
      
    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
    INTEGER(INTG), INTENT(IN) :: datapointusernumber
    INTEGER(INTG), INTENT(IN) :: coupledmeshindexnumber
    INTEGER(INTG), INTENT(OUT) :: coupledmeshuserelementnumber
    INTEGER(INTG), INTENT(IN) :: meshcomponentnumber
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(mesh_type), POINTER :: interfacemesh
    INTEGER(INTG) :: datapointglobalnumber
    LOGICAL :: datapointexists

    enters("InterfacePointsConnectivity_ElementNumberGet",err,error,*999)
    
    IF(ASSOCIATED(pointsconnectivity)) THEN 
      CALL data_point_check_exists(pointsconnectivity%interface%DATA_POINTS,datapointusernumber,datapointexists, &
        & datapointglobalnumber,err,error,*999)
      IF(datapointexists) THEN
        interfacemesh=>pointsconnectivity%interfaceMesh
        coupledmeshuserelementnumber=pointsconnectivity%pointsConnectivity(datapointglobalnumber,coupledmeshindexnumber)% &
          & coupledmeshelementnumber
      ELSE
        CALL flagerror("Data point with user number ("//trim(number_to_vstring &
          & (datapointusernumber,"*",err,error))//") does not exist.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_ElementNumberGet")
    RETURN
999 errors("InterfacePointsConnectivity_ElementNumberGet",err,error)
    exits("InterfacePointsConnectivity_ElementNumberGet")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_elementnumberget
  
  !
  !================================================================================================================================
  !
  
  SUBROUTINE interfacepointsconnectivity_elementnumberset(pointsConnectivity,dataPointUserNumber,coupledMeshIndexNumber, &
      & coupledmeshuserelementnumber,meshcomponentnumber,err,error,*)
      
    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
    INTEGER(INTG), INTENT(IN) :: datapointusernumber
    INTEGER(INTG), INTENT(IN) :: coupledmeshindexnumber
    INTEGER(INTG), INTENT(IN) :: coupledmeshuserelementnumber
    INTEGER(INTG), INTENT(IN) :: meshcomponentnumber
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointglobalnumber,elementmeshnumber
    LOGICAL :: datapointexists,elementexists

    enters("InterfacePointsConnectivity_ElementNumberSet",err,error,*999)
    
    IF(ASSOCIATED(pointsconnectivity)) THEN
      IF(pointsconnectivity%pointsConnectivityFinished) THEN
        CALL flagerror("Interface points connectivity has already been finished.",err,error,*999)
      ELSE
        CALL data_point_check_exists(pointsconnectivity%interface%DATA_POINTS,datapointusernumber,datapointexists, &
          & datapointglobalnumber,err,error,*999)
        IF(datapointexists) THEN
          IF ((coupledmeshindexnumber<=pointsconnectivity%interface%DATA_POINTS%NUMBER_OF_DATA_POINTS).OR. &
              & (coupledmeshindexnumber>0)) THEN
            IF (ALLOCATED(pointsconnectivity%pointsConnectivity)) THEN
              CALL meshtopologyelementcheckexists(pointsconnectivity%INTERFACE%COUPLED_MESHES(coupledmeshindexnumber)%PTR, &
                & meshcomponentnumber,coupledmeshuserelementnumber,elementexists,elementmeshnumber,err,error,*999) !Make sure user element exists       
              IF(elementexists) THEN
                pointsconnectivity%pointsConnectivity(datapointglobalnumber,coupledmeshindexnumber)%coupledMeshElementNumber= &
                  & elementmeshnumber
              ELSE
                CALL flagerror("Element with user number ("//trim(number_to_vstring &
                  & (coupledmeshuserelementnumber,"*",err,error))//") does not exist.",err,error,*999)
              ENDIF
            ELSE
              CALL flagerror("Interface points connectivity array not allocated.",err,error,*999)
            END IF
          ELSE
            CALL flagerror("Interface coupled mesh index number out of range.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Data point with user number ("//trim(number_to_vstring &
              & (datapointusernumber,"*",err,error))//") does not exist.",err,error,*999)
        ENDIF
      ENDIF  
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfacePointsConnectivity_ElementNumberSet")
    RETURN
999 errors("InterfacePointsConnectivity_ElementNumberSet",err,error)
    exits("InterfacePointsConnectivity_ElementNumberSet")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_elementnumberset
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_finalise(interfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: coupledmeshidx,datapointidx
    
    enters("InterfacePointsConnectivity_Finalise",err,error,*999)
    
    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      IF(ALLOCATED(interfacepointsconnectivity%PointsConnectivity)) THEN
        DO coupledmeshidx=1,size(interfacepointsconnectivity%pointsConnectivity,2)
          DO datapointidx=1,size(interfacepointsconnectivity%pointsConnectivity,1) !Deallocate memory for each data point
            CALL interfacepointsconnectivity_pointfinalise(interfacepointsconnectivity%pointsConnectivity(datapointidx, &
              & coupledmeshidx),err,error,*999) 
          ENDDO
        ENDDO
        DEALLOCATE(interfacepointsconnectivity%pointsConnectivity)  
      ENDIF
      CALL interfacepointsconnectivity_coupledelementsfinalise(interfacepointsconnectivity,err,error,*999) 
      NULLIFY(interfacepointsconnectivity%interface)
      NULLIFY(interfacepointsconnectivity%interfaceMesh) 
      IF(ALLOCATED(interfacepointsconnectivity%coupledElements)) DEALLOCATE(interfacepointsconnectivity%coupledElements)
      IF(ALLOCATED(interfacepointsconnectivity%maxNumberOfCoupledElements))  &
        & DEALLOCATE(interfacepointsconnectivity%maxNumberOfCoupledElements)
      DEALLOCATE(interfacepointsconnectivity)
    ENDIF
    
    exits("InterfacePointsConnectivity_Finalise")
    RETURN
999 errorsexits("InterfacePointsConnectivity_Finalise",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_finalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_pointfinalise(interfacePointConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointconnectivitytype) :: interfacepointconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables

    enters("InterfacePointsConnectivity_PointFinalise",err,error,*999)
    
    interfacepointconnectivity%coupledMeshElementNumber=0
    interfacepointconnectivity%elementLineFaceNumber=0
    IF(ALLOCATED(interfacepointconnectivity%xi)) DEALLOCATE(interfacepointconnectivity%xi)
    IF(ALLOCATED(interfacepointconnectivity%reducedXi)) DEALLOCATE(interfacepointconnectivity%reducedXi)
    
    exits("InterfacePointsConnectivity_PointFinalise")
    RETURN
999 errorsexits("InterfacePointsConnectivity_PointFinalise",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_pointfinalise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_fullxicalculate(InterfacePointsConnectivity,coupledMeshIdx,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(IN) :: coupledmeshidx
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointidx,interfacemeshdimensions,coupledmeshdimensions
    TYPE(interface_type), POINTER :: interface
    TYPE(interfacepointconnectivitytype), POINTER :: pointconnectivity
    TYPE(mesh_type), POINTER :: interfacemesh,coupledmesh
    TYPE(varying_string) :: localerror
  
    enters("InterfacePointsConnectivity_FullXiCalculate",err,error,*999)

     IF(ASSOCIATED(interfacepointsconnectivity)) THEN
       interfacemesh=>interfacepointsconnectivity%interfaceMesh
       IF(ASSOCIATED(interfacemesh)) THEN
         interfacemeshdimensions=interfacepointsconnectivity%interfaceMesh%NUMBER_OF_DIMENSIONS
         interface=>interfacepointsconnectivity%interface
         IF(ASSOCIATED(interface)) THEN
           coupledmesh=>interface%COUPLED_MESHES(coupledmeshidx)%PTR
           IF(ASSOCIATED(coupledmesh)) THEN
             coupledmeshdimensions=coupledmesh%NUMBER_OF_DIMENSIONS
             IF(interfacemeshdimensions==coupledmeshdimensions) THEN !e.g. If 1D-2D, 2D-3D coupling, interface dimension is 1D and 2D respectively for 1st body
               DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                 interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)%xi(:)= &
                   & interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)%reducedXi(:)
               ENDDO !dataPointIdx
             ELSE
               !Update full xi location from reduced xi and element face/line number
               SELECT CASE(coupledmeshdimensions)
               CASE(2)
                 DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                   pointconnectivity=>interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)
                   SELECT CASE(pointconnectivity%elementLineFaceNumber)
                   CASE(1)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=0.0_dp
                   CASE(2)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=1.0_dp
                   CASE(3)
                     pointconnectivity%xi(1)=0.0_dp
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(1)
                   CASE(4)
                     pointconnectivity%xi(1)=1.0_dp
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(1)
                   CASE DEFAULT
                     localerror="Invalid local line number "// &
                       & trim(number_to_vstring(pointconnectivity%elementLineFaceNumber,"*",err,error))//" ."
                     CALL flagerror(localerror,err,error,*999)
                   END SELECT
                 ENDDO !dataPointIdx
               CASE(3)
                 DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                   pointconnectivity=>interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshidx)
                   SELECT CASE(pointconnectivity%elementLineFaceNumber)
                   CASE(1)
                     pointconnectivity%xi(1)=1.0_dp
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(3)=pointconnectivity%reducedXi(2)
                   CASE(2)
                     pointconnectivity%xi(1)=0.0_dp
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(3)=pointconnectivity%reducedXi(2)
                   CASE(3)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=1.0_dp
                     pointconnectivity%xi(3)=pointconnectivity%reducedXi(2)
                   CASE(4)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=0.0_dp
                     pointconnectivity%xi(3)=pointconnectivity%reducedXi(2)
                   CASE(5)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(2)
                     pointconnectivity%xi(3)=1.0_dp
                   CASE(6)
                     pointconnectivity%xi(1)=pointconnectivity%reducedXi(1)
                     pointconnectivity%xi(2)=pointconnectivity%reducedXi(2)
                     pointconnectivity%xi(3)=0.0_dp
                   CASE DEFAULT
                     localerror="Invalid local face number "// &
                       & trim(number_to_vstring(pointconnectivity%elementLineFaceNumber,"*",err,error))//" ."
                     CALL flagerror(localerror,err,error,*999)
                   END SELECT
                 ENDDO !dataPointIdx
               CASE DEFAULT
                 localerror="Invalid coupled mesh dimension "// &
                   & trim(number_to_vstring(coupledmeshdimensions,"*",err,error))//" ."
                 CALL flagerror(localerror,err,error,*999)
               END SELECT
             ENDIF
           ELSE
             CALL flagerror("Coupled mesh is not associated.",err,error,*999)
           ENDIF
         ELSE
           CALL flagerror("Interface is not associated.",err,error,*999)
         ENDIF
       ELSE
         CALL flagerror("Interface mesh is not associated.",err,error,*999)
       ENDIF
     ELSE
       CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
     ENDIF
    
    exits("InterfacePointsConnectivity_FullXiCalculate")
    RETURN
999 errorsexits("InterfacePointsConnectivity_FullXiCalculate",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_fullxicalculate
  
    !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_initialise(interface,interfaceMesh,err,error,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    TYPE(mesh_type), POINTER :: interfacemesh
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointidx,meshidx,coupledmeshdimension,interfacemeshdimension
    INTEGER(INTG) :: dummyerr
    TYPE(varying_string)  :: dummyerror
     
    enters("InterfacePointsConnectivity_Initialise",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(ASSOCIATED(interface%pointsConnectivity)) THEN
        CALL flagerror("Interface has already got points connectivity associated.",err,error,*998)
      ELSE
        !Initialise the poins connectivity
        ALLOCATE(interface%pointsConnectivity,stat=err)
        IF(err/=0) CALL flagerror("Could not allocate interface points connectivity.",err,error,*998)
        interface%pointsConnectivity%interface=>interface
        interface%pointsConnectivity%pointsConnectivityFinished=.false.
        interface%pointsConnectivity%interfaceMesh=>interfacemesh
        interfacemeshdimension=interfacemesh%NUMBER_OF_DIMENSIONS
        IF(interface%NUMBER_OF_COUPLED_MESHES>0) THEN 
          IF(ASSOCIATED(interface%DATA_POINTS)) THEN
            IF(interface%DATA_POINTS%NUMBER_OF_DATA_POINTS>0) THEN
              ALLOCATE(interface%pointsConnectivity%pointsConnectivity(interface%DATA_POINTS%NUMBER_OF_DATA_POINTS, &
                & interface%NUMBER_OF_COUPLED_MESHES),stat=err)
              IF(err/=0) CALL flagerror("Could not allocate interface point connectivity.",err,error,*999)
              DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                DO meshidx=1,interface%NUMBER_OF_COUPLED_MESHES
                  coupledmeshdimension=interface%COUPLED_MESHES(meshidx)%PTR%NUMBER_OF_DIMENSIONS
                  CALL interfacepointsconnectivity_pointinitialise(interface%pointsConnectivity%pointsConnectivity(datapointidx, &
                    & meshidx),coupledmeshdimension,interfacemeshdimension,err,error,*999)
                ENDdo!meshIdx
              ENDdo!dataPointIdx
            ELSE
              CALL flagerror("Number of interface data points must be > 0.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Interface data points are not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Number of coupled meshes in the interface must be > 0.",err,error,*999)
        ENDIF 
        ALLOCATE(interface%pointsConnectivity%maxNumberOfCoupledElements(interface%NUMBER_OF_COUPLED_MESHES),stat=err)
        IF(err/=0) CALL flagerror("Could not allocate interface max number of coupled mesh elements.",err,error,*999)
        interface%pointsConnectivity%maxNumberOfCoupledElements=0
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF

    exits("InterfacePointsConnectivity_Initialise")
    RETURN
999 CALL interfacepointsconnectivity_finalise(interface%pointsConnectivity,dummyerr,dummyerror,*998) 
998 errorsexits("InterfacePointsConnectivity_Initialise",err,error)
    RETURN 1
  END SUBROUTINE interfacepointsconnectivity_initialise
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_pointinitialise(interfacePointConnectivity,coupledMeshDimension,interfaceMeshDimension, &
      & err,error,*)

    !Argument variables
    TYPE(interfacepointconnectivitytype) :: interfacepointconnectivity
    INTEGER(INTG), INTENT(IN) :: coupledmeshdimension,interfacemeshdimension
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: dummyerr
    TYPE(varying_string)  :: dummyerror
     
    enters("InterfacePointsConnectivity_PointInitialise",err,error,*999)

    interfacepointconnectivity%coupledMeshElementNumber=0
    interfacepointconnectivity%elementLineFaceNumber=0
    !Allocate memory for coupled mesh full and reduced xi location
    ALLOCATE(interfacepointconnectivity%xi(coupledmeshdimension),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate interface point connectivity full xi.",err,error,*999)
    interfacepointconnectivity%xi=0.0_dp
    ALLOCATE(interfacepointconnectivity%reducedXi(interfacemeshdimension),stat=err)
    IF(err/=0) CALL flagerror("Could not allocate interface point connectivity reduced xi.",err,error,*999)
    interfacepointconnectivity%reducedXi=0.0_dp
       
    exits("InterfacePointsConnectivity_PointInitialise")
    RETURN
999 CALL interfacepointsconnectivity_pointfinalise(interfacepointconnectivity,dummyerr,dummyerror,*998)
998 errorsexits("InterfacePointsConnectivity_PointInitialise",err,error)
    RETURN 1
  END SUBROUTINE interfacepointsconnectivity_pointinitialise
  
  !
  !================================================================================================================================
  !
    
  SUBROUTINE interfacepointsconnectivity_pointxiget(pointsConnectivity,dataPointUserNumber,coupledMeshIndexNumber, &
    & xi,err,error,*)

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
    INTEGER(INTG), INTENT(IN) :: datapointusernumber
    INTEGER(INTG), INTENT(IN) :: coupledmeshindexnumber
    REAL(DP), INTENT(OUT) :: xi(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointglobalnumber
    LOGICAL :: datapointexists
    
    enters("InterfacePointsConnectivity_PointXiGet",err,error,*999)
    
    ! Preliminary error checks to verify user input information
    IF(ASSOCIATED(pointsconnectivity)) THEN
      IF (ALLOCATED(pointsconnectivity%pointsConnectivity)) THEN
        CALL data_point_check_exists(pointsconnectivity%INTERFACE%DATA_POINTS,datapointusernumber,datapointexists, &
          & datapointglobalnumber,err,error,*999)
        IF(datapointexists) THEN
          IF(SIZE(xi)>=SIZE(pointsconnectivity%pointsConnectivity(datapointusernumber,coupledmeshindexnumber)%xi,1)) THEN
            xi=pointsconnectivity%pointsConnectivity(datapointusernumber,coupledmeshindexnumber)%xi(:)
          ELSE
            CALL flagerror("Input xi array size is smaller than points connectivity xi array.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Data point with user number ("//trim(number_to_vstring &
              & (datapointusernumber,"*",err,error))//") does not exist.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface points connectivity array not allocated.",err,error,*999)
      ENDIF 
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF   

    exits("InterfacePointsConnectivity_PointXiGet")
    RETURN
999 errorsexits("InterfacePointsConnectivity_PointXiGet",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_pointxiget
  
  !
  !================================================================================================================================
  !
    
  SUBROUTINE interfacepointsconnectivity_pointxiset(pointsConnectivity,dataPointUserNumber,coupledMeshIndexNumber, &
    & xi,err,error,*)

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: pointsconnectivity
    INTEGER(INTG), INTENT(IN) :: datapointusernumber
    INTEGER(INTG), INTENT(IN) :: coupledmeshindexnumber
    REAL(DP), INTENT(IN) :: xi(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointglobalnumber
    LOGICAL :: datapointexists
    
    enters("InterfacePointsConnectivity_PointXiSet",err,error,*999)
    
    ! Preliminary error checks to verify user input information
    IF(ASSOCIATED(pointsconnectivity)) THEN
      IF(pointsconnectivity%pointsConnectivityFinished) THEN
        CALL flagerror("Interface mesh connectivity already been finished.",err,error,*999)
      ELSE
        IF (ALLOCATED(pointsconnectivity%pointsConnectivity)) THEN
          CALL data_point_check_exists(pointsconnectivity%INTERFACE%DATA_POINTS,datapointusernumber,datapointexists, &
            & datapointglobalnumber,err,error,*999)
          IF(datapointexists) THEN
            IF ((coupledmeshindexnumber>pointsconnectivity%interface%NUMBER_OF_COUPLED_MESHES).OR.(coupledmeshindexnumber<0)) THEN
              CALL flagerror("Interface coupled mesh index number out of range.",err,error,*999)
            ELSE
              IF(SIZE(pointsconnectivity%pointsConnectivity(datapointglobalnumber,coupledmeshindexnumber)%xi,1)== &
                  & SIZE(xi,1)) THEN
                pointsconnectivity%pointsConnectivity(datapointglobalnumber,coupledmeshindexnumber)% &
                  & xi(:)=xi(:) 
              ELSE
                CALL flagerror("Input xi dimension does not match full coupled mesh xi dimension.",err,error,*999)
              ENDIF
            ENDIF
          ELSE
            CALL flagerror("Data point with user number ("//trim(number_to_vstring &
                & (datapointusernumber,"*",err,error))//") does not exist.",err,error,*999)
          ENDIF
        ELSE
          CALL flagerror("Interface points connectivity array not allocated.",err,error,*999)
        ENDIF 
      ENDIF
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF   

    exits("InterfacePointsConnectivity_PointXiSet")
    RETURN
999 errorsexits("InterfacePointsConnectivity_PointXiSet",err,error)
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_pointxiset
  
  !
  !================================================================================================================================
  !
  
  SUBROUTINE interfacepointsconnectivity_updatefromprojection(InterfacePointsConnectivity,dataProjection, &
      & coupledmeshindex,err,error,*) 
  
    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    TYPE(data_projection_type), POINTER :: dataprojection
    INTEGER(INTG), INTENT(IN) :: coupledmeshindex
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: datapointidx
    TYPE(data_projection_result_type), POINTER :: dataprojectionresult
    
    enters("InterfacePointsConnectivity_UpdateFromProjection",err,error,*999)
    
    IF(ASSOCIATED(interfacepointsconnectivity)) THEN
      IF(ASSOCIATED(dataprojection)) THEN
        IF(dataprojection%DATA_PROJECTION_FINISHED) THEN
          WRITE(*,*) "InterfacePointsConnectivity_UpdateFromProjection"
          DO datapointidx=1,SIZE(dataprojection%DATA_PROJECTION_RESULTS,1) !Update reduced xi location, projection element number and element face/line number with projection result
            dataprojectionresult=>dataprojection%DATA_PROJECTION_RESULTS(datapointidx)
            interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshindex)%reducedXi(:)= &
              & dataprojectionresult%XI
            interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshindex)%coupledMeshElementNumber= &
              & dataprojectionresult%ELEMENT_NUMBER
            IF(dataprojectionresult%ELEMENT_LINE_NUMBER/=0) THEN
              interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshindex)%elementLineFaceNumber= &
                & dataprojectionresult%ELEMENT_LINE_NUMBER
            ELSEIF(dataprojectionresult%ELEMENT_FACE_NUMBER/=0) THEN
              interfacepointsconnectivity%pointsConnectivity(datapointidx,coupledmeshindex)%elementLineFaceNumber= &
                & dataprojectionresult%ELEMENT_FACE_NUMBER
            ENDIF
          ENDDO
          CALL interfacepointsconnectivity_fullxicalculate(interfacepointsconnectivity,coupledmeshindex, &
            & err,error,*999) 
          !Update points connectivity coupledElement information
          CALL interfacepointsconnectivity_coupledelementscalculate(interfacepointsconnectivity,coupledmeshindex,err,error,*999) 
        ELSE
          CALL flagerror("Data projection is not finished.",err,error,*999)
        ENDIF
      ELSE
        CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
    ENDIF
     
    exits("InterfacePointsConnectivity_UpdateFromProjection")
    RETURN
999 errors("InterfacePointsConnectivity_UpdateFromProjection",err,error)
    exits("InterfacePointsConnectivity_UpdateFromProjection")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_updatefromprojection
  
  !
  !================================================================================================================================
  !

  SUBROUTINE interfacepointsconnectivity_reducedxicalculate(InterfacePointsConnectivity,err,error,*) 

    !Argument variables
    TYPE(interfacepointsconnectivitytype), POINTER :: interfacepointsconnectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: meshidx,datapointidx,xiidx,interfacemeshdimensions,coupledmeshdimensions
    TYPE(interface_type), POINTER :: interface
  
    enters("InterfacePointsConnectivity_ReducedXiCalculate",err,error,*999)

     IF(ASSOCIATED(interfacepointsconnectivity)) THEN
       IF(interfacepointsconnectivity%pointsConnectivityFinished) THEN
         CALL flagerror("Interface points connectivity has already been finished.",err,error,*999)
       ELSE
         interface=>interfacepointsconnectivity%interface
         interfacemeshdimensions=interfacepointsconnectivity%interfaceMesh%NUMBER_OF_DIMENSIONS
         IF(ASSOCIATED(interface)) THEN
           DO meshidx=1,interface%NUMBER_OF_COUPLED_MESHES
             coupledmeshdimensions=interface%COUPLED_MESHES(meshidx)%PTR%NUMBER_OF_DIMENSIONS
             IF(interfacemeshdimensions==coupledmeshdimensions) THEN !e.g. If 1D-2D, 2D-3D coupling, interface dimension is 1D and 2D respectively for 1st body
               DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                 interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(:)= &
                   & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(:)
                 interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%elementLineFaceNumber=1 !The only local face/line for the body with lower dimension 
               ENDDO !dataPointIdx
             ELSE
               SELECT CASE(coupledmeshdimensions)
               CASE(2)
                 DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                   DO xiidx=1,coupledmeshdimensions
                     IF(abs(interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(xiidx)) &
                       & < zero_tolerance) THEN
                       interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%elementLineFaceNumber=4-(xiidx-1)*2 !Calculate line number
                       EXIT
                     ELSEIF(abs(interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(xiidx)-1.0_dp) &
                       & < zero_tolerance) THEN
                       interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%elementLineFaceNumber=3-(xiidx-1)*2 !Calculate line number
                       EXIT
                     ENDIF
                   ENDDO !xiIdx
                   interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi= &
                     & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)% &
                     & xi(other_xi_directions2(xiidx))  !Populate reducedXi 
                 ENDDO !dataPointIdx
               CASE(3)
                 DO datapointidx=1,interface%DATA_POINTS%NUMBER_OF_DATA_POINTS
                   DO xiidx=1,coupledmeshdimensions
                     IF(abs(interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(xiidx)) &
                       & < zero_tolerance) THEN
                       interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%elementLineFaceNumber=(xiidx-1)*2+2 !Calculate face number
                       EXIT
                     ELSE IF(abs(interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(xiidx)-1.0_dp) &
                       & > zero_tolerance) THEN
                       interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%elementLineFaceNumber=(xiidx-1)*2+1 !Calculate face number
                       EXIT
                     ENDIF
                   ENDDO !xiIdx
                   SELECT CASE(xiidx) !Populate reducedXi 
                   CASE(1)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(1)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(2)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(2)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(3)
                   CASE(2)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(1)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(1)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(2)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(3)
                   CASE(3)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(1)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(1)
                     interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%reducedXi(2)= &
                      & interfacepointsconnectivity%pointsConnectivity(datapointidx,meshidx)%xi(2)
                   END SELECT
                 ENDDO !dataPointIdx
               CASE DEFAULT
                 ! Do nothing
               END SELECT
             ENDIF
           ENDDO !meshIdx
         ELSE
           CALL flagerror("Interface is not associated.",err,error,*999)
         ENDIF
       ENDIF
     ELSE
       CALL flagerror("Interface points connectivity is not associated.",err,error,*999)
     ENDIF
    
    exits("InterfacePointsConnectivity_ReducedXiCalculate")
    RETURN
999 errors("InterfacePointsConnectivity_ReducedXiCalculate",err,error)
    exits("InterfacePointsConnectivity_ReducedXiCalculate")
    RETURN 1
    
  END SUBROUTINE interfacepointsconnectivity_reducedxicalculate

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

  SUBROUTINE interface_user_number_find(USER_NUMBER,PARENT_REGION,INTERFACE,ERR,ERROR,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: user_number
    TYPE(region_type), POINTER :: parent_region
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interface_idx
    TYPE(varying_string) :: local_error

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

    IF(ASSOCIATED(parent_region)) THEN
      IF(ASSOCIATED(interface)) THEN
        CALL flagerror("Interface is already associated.",err,error,*999)
      ELSE
        NULLIFY(interface)
        IF(ASSOCIATED(parent_region%INTERFACES)) THEN
          interface_idx=1
          DO WHILE(interface_idx<=parent_region%INTERFACES%NUMBER_OF_INTERFACES.AND..NOT.ASSOCIATED(interface))
            IF(parent_region%INTERFACES%INTERFACES(interface_idx)%PTR%USER_NUMBER==user_number) THEN
              interface=>parent_region%INTERFACES%INTERFACES(interface_idx)%PTR
            ELSE
              interface_idx=interface_idx+1
            ENDIF
          ENDDO
        ELSE
          local_error="The interfaces on parent region number "// &
            & trim(number_to_vstring(parent_region%USER_NUMBER,"*",err,error))//" are not associated."
          CALL flagerror(local_error,err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flagerror("Parent region is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACE_USER_NUMBER_FIND")
    RETURN
999 errorsexits("INTERFACE_USER_NUMBER_FIND",err,error)
    RETURN 1
  END SUBROUTINE interface_user_number_find

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

  SUBROUTINE interfaces_finalise(INTERFACES,ERR,ERROR,*) 

    !Argument variables
    TYPE(interfaces_type), POINTER :: interfaces
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(interface_type), POINTER :: interface
     
    enters("INTERFACES_FINALISE",err,error,*999)

    IF(ASSOCIATED(interfaces)) THEN
      DO WHILE(interfaces%NUMBER_OF_INTERFACES>0)
        interface=>interfaces%INTERFACES(1)%PTR
        CALL interface_destroy(interface,err,error,*999)
      ENDDO
      IF(ASSOCIATED(interfaces%INTERFACES)) DEALLOCATE(interfaces%INTERFACES)
      DEALLOCATE(interfaces)
    ENDIF
    
    exits("INTERFACES_FINALISE")
    RETURN
999 errorsexits("INTERFACES_FINALISE",err,error)
    RETURN 1
  END SUBROUTINE interfaces_finalise

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

  SUBROUTINE interfaces_initialise(REGION,ERR,ERROR,*) 

    !Argument variables
    TYPE(region_type), POINTER :: region
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(varying_string) :: local_error
     
    enters("INTERFACES_INITIALISE",err,error,*998)

    IF(ASSOCIATED(region)) THEN
      IF(ASSOCIATED(region%INTERFACES)) THEN
        local_error="Region number "//trim(number_to_vstring(region%USER_NUMBER,"*",err,error))// &
          & " already has interfaces associated."
        CALL flagerror(local_error,err,error,*999)
      ELSE
        ALLOCATE(region%INTERFACES,stat=err)
        IF(err/=0) CALL flagerror("Could not allocate region interfaces.",err,error,*999)
        region%INTERFACES%PARENT_REGION=>region
        region%INTERFACES%NUMBER_OF_INTERFACES=0
        NULLIFY(region%INTERFACES%INTERFACES)
      ENDIF
    ELSE
      CALL flagerror("Region is not associated.",err,error,*999)
    ENDIF
    
    exits("INTERFACES_INITIALISE")
    RETURN
999 CALL interfaces_finalise(region%INTERFACES,err,error,*998)
998 errorsexits("INTERFACES_INITIALISE",err,error)
    RETURN 1
  END SUBROUTINE interfaces_initialise

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

  SUBROUTINE interfacemeshconnectivity_elementinitialise(INTERFACE,ERR,ERROR,*)

    !Argument variables
    TYPE(interface_type), POINTER :: interface
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interfaceelementidx,coupledmeshidx
    TYPE(interface_element_connectivity_type), POINTER :: element_connectivity
     
    enters("InterfaceMeshConnectivity_ElementInitialise",err,error,*999)

    IF(ASSOCIATED(interface)) THEN
      IF(ASSOCIATED(interface%MESH_CONNECTIVITY)) THEN
        IF(ALLOCATED(interface%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY)) THEN
          CALL flagerror("Interface mesh element connectivity is already allocated.",err,error,*999)
        ELSE
          IF(interface%NUMBER_OF_COUPLED_MESHES>0) THEN
            IF(interface%MESH_CONNECTIVITY%INTERFACE_MESH%NUMBER_OF_ELEMENTS>0) THEN
              interface%MESH_CONNECTIVITY%NUMBER_OF_INTERFACE_ELEMENTS=interface%MESHES%MESHES(1)%PTR%NUMBER_OF_ELEMENTS
              interface%MESH_CONNECTIVITY%NUMBER_OF_COUPLED_MESHES=interface%NUMBER_OF_COUPLED_MESHES
              ALLOCATE(interface%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(interface%MESH_CONNECTIVITY%NUMBER_OF_INTERFACE_ELEMENTS, &
                & interface%MESH_CONNECTIVITY%NUMBER_OF_COUPLED_MESHES),stat=err)
              IF(err/=0) CALL flagerror("Could not allocate interface element connectivity.",err,error,*999)
              DO interfaceelementidx=1,interface%MESH_CONNECTIVITY%NUMBER_OF_INTERFACE_ELEMENTS
                DO coupledmeshidx=1,interface%MESH_CONNECTIVITY%NUMBER_OF_COUPLED_MESHES
                  element_connectivity=>interface%MESH_CONNECTIVITY%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)
                  element_connectivity%COUPLED_MESH_ELEMENT_NUMBER=0
                  element_connectivity%CONNECTED_FACE=0
                  element_connectivity%CONNECTED_LINE=0
                  ! Note that ELEMENT_CONNECTIVITY%XI(NumberOfCoupledMeshXiDirections,1,NumberOfInterfaceElementNodes) is allocated after the basis for the mesh connectivity has been specified in INTERFACE_MESH_CONNECTIVITY_BASIS_SET where the number of NumberOfInterfaceElementNodes can be determined.
                  !\todo see corresponding todo in regarding updating the structure of ELEMENT_CONNECTIVITY%XI
                ENDDO !CoupledMeshIdx
              ENDDO !InterfaceElementIdx
            ELSE
              CALL flagerror("Interface coupled meshes are not associated.",err,error,*999)
            ENDIF
          ELSE
            CALL flagerror("Interface coupled meshes are not associated.",err,error,*999)
          ENDIF
        ENDIF
      ELSE
        CALL flagerror("Interface mesh connectivity is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flagerror("Interface is not associated.",err,error,*999)
    ENDIF
    
    exits("InterfaceMeshConnectivity_ElementInitialise")
    RETURN
999 errorsexits("InterfaceMeshConnectivity_ElementInitialise",err,error)
    RETURN 1
    
  END SUBROUTINE interfacemeshconnectivity_elementinitialise

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

  SUBROUTINE interfacemeshconnectivity_elementfinalise(INTERFACE_MESH_CONNECTIVITY,ERR,ERROR,*) 

    !Argument variables
    TYPE(interface_mesh_connectivity_type) :: interface_mesh_connectivity
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: interfaceelementidx,coupledmeshidx
     
    enters("InterfaceMeshConnectivity_ElementFinalise",err,error,*999)

    DO interfaceelementidx=1,interface_mesh_connectivity%NUMBER_OF_INTERFACE_ELEMENTS  
      DO coupledmeshidx=1,interface_mesh_connectivity%NUMBER_OF_COUPLED_MESHES
        IF(ALLOCATED(interface_mesh_connectivity%ELEMENT_CONNECTIVITY)) THEN
          interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)%COUPLED_MESH_ELEMENT_NUMBER=0
          IF(ALLOCATED(interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)%XI)) THEN
            DEALLOCATE(interface_mesh_connectivity%ELEMENT_CONNECTIVITY(interfaceelementidx,coupledmeshidx)%XI)
          ENDIF
        ELSE
          CALL flagerror("Interface mesh connectivity element connectivity is being deallocated before allocation.", &
            & err,error,*999)
        ENDIF
      ENDDO !InterfaceElementIdx
    ENDDO !CoupledMeshIdx

    DEALLOCATE(interface_mesh_connectivity%ELEMENT_CONNECTIVITY)
    
    exits("InterfaceMeshConnectivity_ElementFinalise")
    RETURN
999 errorsexits("InterfaceMeshConnectivity_ElementFinalise",err,error)
    RETURN 1
    
  END SUBROUTINE interfacemeshconnectivity_elementfinalise

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

END MODULE interface_routines