data_projection_routines.f90

Go to the documentation of this file.

MODULE data_projection_routines

  USE base_routines
  USE cmiss_mpi
  USE comp_environment
  USE constants
  USE domain_mappings
  USE field_routines
  USE input_output
  USE iso_varying_string
  USE kinds
  USE mesh_routines
#ifndef NOMPIMOD
  USE mpi
#endif
  USE sorting
  USE strings
  USE trees
  USE types

#include "macros.h"

  IMPLICIT NONE

  PRIVATE

#ifdef NOMPIMOD
#include "mpif.h"
#endif

  !Module parameters

  INTEGER(INTG), PARAMETER :: data_projection_boundary_lines_projection_type=1
  INTEGER(INTG), PARAMETER :: data_projection_boundary_faces_projection_type=2
  INTEGER(INTG), PARAMETER :: data_projection_all_elements_projection_type=3

  INTEGER(INTG), PARAMETER :: data_projection_exit_tag_converged=1
  INTEGER(INTG), PARAMETER :: data_projection_exit_tag_bounds=2
  INTEGER(INTG), PARAMETER :: data_projection_exit_tag_max_iteration=3
  INTEGER(INTG), PARAMETER :: data_projection_exit_tag_no_element=4

  !Module types

  !Module variables

  !Interfaces

  INTERFACE data_projection_label_get
    MODULE PROCEDURE data_projection_label_get_c
    MODULE PROCEDURE data_projection_label_get_vs
  END INTERFACE !DATA_PROJECTION_LABEL_GET

  INTERFACE data_projection_label_set
    MODULE PROCEDURE data_projection_label_set_c
    MODULE PROCEDURE data_projection_label_set_vs
  END INTERFACE !DATA_PROJECTION_LABEL_SET

  PUBLIC data_projection_boundary_lines_projection_type,data_projection_boundary_faces_projection_type, &
    & data_projection_all_elements_projection_type

  PUBLIC data_projection_absolute_tolerance_get,data_projection_absolute_tolerance_set

  PUBLIC data_projection_create_finish,data_projection_create_start_data_points

  PUBLIC data_projection_destroy

  PUBLIC dataprojection_datapointsprojectionevaluate

  PUBLIC dataprojection_datapointspositionevaluate

  PUBLIC dataprojection_maximuminterationupdateget,dataprojection_maximuminterationupdateset

  PUBLIC dataprojection_maximumnumberofiterationsget,dataprojection_maximumnumberofiterationsset

  PUBLIC dataprojection_numberofclosestelementsget,dataprojection_numberofclosestelementsset

  PUBLIC dataprojection_projectioncandidatesset

  PUBLIC data_projection_projection_type_get,data_projection_projection_type_set

  PUBLIC data_projection_relative_tolerance_get,data_projection_relative_tolerance_set

  PUBLIC data_projection_starting_xi_get,data_projection_starting_xi_set

  PUBLIC data_projection_result_xi_get, data_projection_result_xi_set

  PUBLIC dataprojection_resultprojectionvectorget

  PUBLIC data_projection_element_set

  PUBLIC data_projection_result_distance_get

  PUBLIC data_projection_result_element_number_get,dataprojection_resultelementfacenumberget

  PUBLIC dataprojection_resultelementlinenumberget

  PUBLIC data_projection_result_exit_tag_get

  PUBLIC data_projection_label_get,data_projection_label_set

CONTAINS

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

  SUBROUTINE data_projection_absolute_tolerance_get(DATA_PROJECTION,ABSOLUTE_TOLERANCE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(OUT) :: ABSOLUTE_TOLERANCE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_ABSOLUTE_TOLERANCE_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_ABSOLUTE_TOLERANCE_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_absolute_tolerance_get

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

  SUBROUTINE data_projection_absolute_tolerance_set(DATA_PROJECTION,ABSOLUTE_TOLERANCE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(IN) :: ABSOLUTE_TOLERANCE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF(absolute_tolerance>=0) THEN
          data_projection%ABSOLUTE_TOLERANCE=absolute_tolerance
        ELSE
          CALL flag_error("Data projection absolute tolerance must be a positive real number.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_ABSOLUTE_TOLERANCE_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_ABSOLUTE_TOLERANCE_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_absolute_tolerance_set

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

  SUBROUTINE data_projection_closest_elements_find(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES, &
    & candidate_elements,number_of_candidates,closest_elements,closest_distances,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_CANDIDATES
    INTEGER(INTG), INTENT(OUT) :: CLOSEST_ELEMENTS(:)
    REAL(DP), INTENT(OUT) :: CLOSEST_DISTANCES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

    INTEGER(INTG) :: REGION_DIMENSIONS
    INTEGER(INTG) :: NUMBER_OF_CLOSEST_CANDIDATES
    REAL(DP) :: DISTANCE_VECTOR(3)
    REAL(DP) :: DISTANCE2
    INTEGER(INTG) :: nce
    INTEGER(INTG) :: ELEMENT_NUMBER,insert_idx

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        number_of_closest_candidates=SIZE(closest_elements,1)
        !loop through the first few elements
        DO nce=1,number_of_closest_candidates
          element_number=candidate_elements(nce)
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions) = point_values-interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          DO insert_idx=nce,1,-1
            IF(insert_idx>1) THEN
              IF(distance2<closest_distances(insert_idx-1)) cycle
            ENDIF
            !insert the element number into the correct index
            IF(insert_idx<nce) THEN
              closest_distances(insert_idx+1:nce)=closest_distances(insert_idx:nce-1)
              closest_elements(insert_idx+1:nce)=closest_elements(insert_idx:nce-1)
            ENDIF
            closest_distances(insert_idx)=distance2
            closest_elements(insert_idx)=element_number
            EXIT
          ENDDO
        ENDDO
        !Loop through the rest of the elements
        DO nce=number_of_closest_candidates+1,number_of_candidates
          element_number=candidate_elements(nce)
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values - interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          IF(distance2<closest_distances(number_of_closest_candidates))THEN
            DO insert_idx=number_of_closest_candidates,1,-1
              IF(insert_idx>1) THEN
                IF(distance2<closest_distances(insert_idx-1)) cycle
              ENDIF
              !insert the element into the correct index
              IF(insert_idx<number_of_closest_candidates) THEN
                closest_distances(insert_idx+1:number_of_closest_candidates)=closest_distances(insert_idx: &
                  & number_of_closest_candidates-1)
                closest_elements(insert_idx+1:number_of_closest_candidates)=closest_elements(insert_idx: &
                  & number_of_closest_candidates-1)
              ENDIF
              closest_distances(insert_idx)=distance2
              closest_elements(insert_idx)=element_number
              EXIT
            ENDDO
          ENDIF
        ENDDO
        !CLOSEST_DISTANCES=SQRT(CLOSEST_DISTANCES) !return shortest distances
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_CLOSEST_ELEMENTS_FIND")
    RETURN
999 errorsexits("DATA_PROJECTION_CLOSEST_ELEMENTS_FIND",err,error)
    RETURN 1

  END SUBROUTINE data_projection_closest_elements_find

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

  SUBROUTINE data_projection_closest_faces_find(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & candidate_element_faces,number_of_candidates,closest_elements,closest_element_faces,closest_distances,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENT_FACES(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_CANDIDATES
    INTEGER(INTG), INTENT(OUT) :: CLOSEST_ELEMENTS(:)
    INTEGER(INTG), INTENT(OUT) :: CLOSEST_ELEMENT_FACES(:)
    REAL(DP), INTENT(OUT) :: CLOSEST_DISTANCES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: REGION_DIMENSIONS
    INTEGER(INTG) :: NUMBER_OF_CLOSEST_CANDIDATES
    REAL(DP) :: DISTANCE_VECTOR(3)
    REAL(DP) :: DISTANCE2
    INTEGER(INTG) :: nce
    INTEGER(INTG) :: ELEMENT_NUMBER,ELEMENT_FACE_NUMBER,FACE_NUMBER,insert_idx

    enters("DATA_PROJECTION_CLOSEST_FACES_FIND",err,error,*999)
    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        number_of_closest_candidates=SIZE(closest_elements,1)
        !loop through the first few faces
        DO nce=1,number_of_closest_candidates
          element_number=candidate_elements(nce)
          element_face_number=candidate_element_faces(nce)
          face_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_FACES(element_face_number)
          CALL field_interpolation_parameters_face_get(field_values_set_type,face_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions) = point_values-interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          DO insert_idx=nce,1,-1
            IF(insert_idx>1) THEN
              IF(distance2<closest_distances(insert_idx-1)) cycle
            ENDIF
            !insert the element number into the correct index
            IF(insert_idx<nce) THEN
              closest_distances(insert_idx+1:nce)=closest_distances(insert_idx:nce-1)
              closest_elements(insert_idx+1:nce)=closest_elements(insert_idx:nce-1)
              closest_element_faces(insert_idx+1:nce)=closest_element_faces(insert_idx:nce-1)
            ENDIF
            closest_distances(insert_idx)=distance2
            closest_elements(insert_idx)=element_number
            closest_element_faces(insert_idx)=element_face_number
            EXIT
          ENDDO
        ENDDO
        !Loop through the rest of the faces
        DO nce=number_of_closest_candidates+1,number_of_candidates
          element_number=candidate_elements(nce)
          element_face_number=candidate_element_faces(nce)
          face_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_FACES(element_face_number)
          CALL field_interpolation_parameters_face_get(field_values_set_type,face_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values - interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          IF(distance2<closest_distances(number_of_closest_candidates))THEN
            DO insert_idx=number_of_closest_candidates,1,-1
              IF(insert_idx>1) THEN
                IF(distance2<closest_distances(insert_idx-1)) cycle
              ENDIF
              !insert the element into the correct index
              IF(insert_idx<number_of_closest_candidates) THEN
                closest_distances(insert_idx+1:number_of_closest_candidates)=closest_distances(insert_idx: &
                  & number_of_closest_candidates-1)
                closest_elements(insert_idx+1:number_of_closest_candidates)=closest_elements(insert_idx: &
                  & number_of_closest_candidates-1)
                closest_element_faces(insert_idx+1:number_of_closest_candidates)=closest_element_faces(insert_idx: &
                  & number_of_closest_candidates-1)
              ENDIF
              closest_distances(insert_idx)=distance2
              closest_elements(insert_idx)=element_number
              closest_element_faces(insert_idx)=element_face_number
              EXIT
            ENDDO
          ENDIF
        ENDDO
        !CLOSEST_DISTANCES=SQRT(CLOSEST_DISTANCES) !return shortest distances
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_CLOSEST_FACES_FIND")
    RETURN
999 errorsexits("DATA_PROJECTION_CLOSEST_FACES_FIND",err,error)
    RETURN 1

  END SUBROUTINE data_projection_closest_faces_find

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

  SUBROUTINE data_projection_closest_lines_find(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & candidate_element_lines,number_of_candidates,closest_elements,closest_element_lines,closest_distances,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENT_LINES(:)
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_CANDIDATES
    INTEGER(INTG), INTENT(OUT) :: CLOSEST_ELEMENTS(:)
    INTEGER(INTG), INTENT(OUT) :: CLOSEST_ELEMENT_LINES(:)
    REAL(DP), INTENT(OUT) :: CLOSEST_DISTANCES(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: REGION_DIMENSIONS
    INTEGER(INTG) :: NUMBER_OF_CLOSEST_CANDIDATES
    REAL(DP) :: DISTANCE_VECTOR(3)
    REAL(DP) :: DISTANCE2
    INTEGER(INTG) :: nce
    INTEGER(INTG) :: ELEMENT_NUMBER,ELEMENT_LINE_NUMBER,LINE_NUMBER,insert_idx

    enters("DATA_PROJECTION_CLOSEST_LINES_FIND",err,error,*999)
    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        number_of_closest_candidates=SIZE(closest_elements,1)
        !loop through the first few lines
        DO nce=1,number_of_closest_candidates
          element_number=candidate_elements(nce)
          element_line_number=candidate_element_lines(nce)
          line_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_LINES(element_line_number)
          CALL field_interpolation_parameters_line_get(field_values_set_type,line_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions) = point_values-interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          DO insert_idx=nce,1,-1
            IF(insert_idx>1) THEN
              IF(distance2<closest_distances(insert_idx-1)) cycle
            ENDIF
            !insert the element number into the correct index
            IF(insert_idx<nce) THEN
              closest_distances(insert_idx+1:nce)=closest_distances(insert_idx:nce-1)
              closest_elements(insert_idx+1:nce)=closest_elements(insert_idx:nce-1)
              closest_element_lines(insert_idx+1:nce)=closest_element_lines(insert_idx:nce-1)
            ENDIF
            closest_distances(insert_idx)=distance2
            closest_elements(insert_idx)=element_number
            closest_element_lines(insert_idx)=element_line_number
            EXIT
          ENDDO
        ENDDO
        !Loop through the rest of the lines
        DO nce=number_of_closest_candidates+1,number_of_candidates
          element_number=candidate_elements(nce)
          element_line_number=candidate_element_lines(nce)
          line_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_LINES(element_line_number)
          CALL field_interpolation_parameters_line_get(field_values_set_type,line_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          CALL field_interpolate_xi(no_part_deriv,data_projection%STARTING_XI,interpolated_point,err,error,*999, &
            & field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values - interpolated_point%VALUES(:,1)
          distance2 = dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          IF(distance2<closest_distances(number_of_closest_candidates))THEN
            DO insert_idx=number_of_closest_candidates,1,-1
              IF(insert_idx>1) THEN
                IF(distance2<closest_distances(insert_idx-1)) cycle
              ENDIF
              !insert the element into the correct index
              IF(insert_idx<number_of_closest_candidates) THEN
                closest_distances(insert_idx+1:number_of_closest_candidates)=closest_distances(insert_idx: &
                  & number_of_closest_candidates-1)
                closest_elements(insert_idx+1:number_of_closest_candidates)=closest_elements(insert_idx: &
                  & number_of_closest_candidates-1)
                closest_element_lines(insert_idx+1:number_of_closest_candidates)=closest_element_lines(insert_idx: &
                  & number_of_closest_candidates-1)
              ENDIF
              closest_distances(insert_idx)=distance2
              closest_elements(insert_idx)=element_number
              closest_element_lines(insert_idx)=element_line_number
              EXIT
            ENDDO
          ENDIF
        ENDDO
        !CLOSEST_DISTANCES=SQRT(CLOSEST_DISTANCES) !return shortest distances
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_CLOSEST_LINES_FIND")
    RETURN
999 errorsexits("DATA_PROJECTION_CLOSEST_LINES_FIND",err,error)
    RETURN 1

  END SUBROUTINE data_projection_closest_lines_find

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

  SUBROUTINE data_projection_create_finish(DATA_PROJECTION,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN !Has to be associated
      IF(ASSOCIATED(data_projection%DATA_POINTS)) THEN !Has to be associated
        IF(data_projection%DATA_POINTS%DATA_POINTS_FINISHED) THEN !Has to be finished
          IF(ASSOCIATED(data_projection%MESH)) THEN !Has to be associated
            IF(data_projection%DATA_PROJECTION_FINISHED) THEN !Cannot be finished
              CALL flag_error("Data projection have already been finished.",err,error,*999)
            ELSE
              data_projection%DATA_PROJECTION_FINISHED=.true.
              CALL dataprojection_dataprojectionresultinitialise(data_projection,err,error, &
                & *999) 
            ENDIF
          ELSE
            CALL flag_error("Data projection mesh is not associated.",err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection data points have not been finished.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection data points is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_CREATE_FINISH")
    RETURN
999 errorsexits("DATA_PROJECTION_CREATE_FINISH",err,error)
    RETURN 1

  END SUBROUTINE data_projection_create_finish

  !
  !================================================================================================================================
  !
  SUBROUTINE data_projection_create_start_data_points(DATA_PROJECTION_USER_NUMBER,DATA_POINTS,MESH,DATA_PROJECTION, &
    & err,error,*)

    !Argument variables
    INTEGER(INTG), INTENT(IN) :: DATA_PROJECTION_USER_NUMBER
    TYPE(data_points_type), POINTER :: DATA_POINTS
    TYPE(mesh_type), POINTER :: MESH
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(data_projection_ptr_type), ALLOCATABLE :: NEW_DATA_PROJECTIONS_PTR(:)
    INTEGER(INTG) :: DATA_POINTS_REGION_DIMENSIONS,MESH_REGION_DIMENSIONS,INSERT_STATUS
    INTEGER(INTG) :: xi_idx,data_projection_idx

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

    NULLIFY(data_projection)
    IF(ASSOCIATED(data_points)) THEN !Has to be associated
      IF(data_points%DATA_POINTS_FINISHED) THEN !Has to be finished
        IF(ASSOCIATED(mesh)) THEN !Has to be associated
          IF(ASSOCIATED(data_projection)) THEN !Cannot be associated
            CALL flag_error("Data projection is already associated.",err,error,*999)
          ELSE
            IF(ASSOCIATED(data_points%REGION)) THEN
              data_points_region_dimensions=data_points%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
            ELSEIF(ASSOCIATED(data_points%INTERFACE)) THEN
              data_points_region_dimensions=data_points%INTERFACE%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
            ELSE
              CALL flag_error("Data points is not associated with a region or a interface.",err,error,*999)
            ENDIF
            IF(ASSOCIATED(mesh%REGION)) THEN
              mesh_region_dimensions=mesh%REGION%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
            ELSEIF(ASSOCIATED(mesh%INTERFACE)) THEN
              mesh_region_dimensions=mesh%INTERFACE%COORDINATE_SYSTEM%NUMBER_OF_DIMENSIONS
            ELSE
              CALL flag_error("Mesh is not associated with a region or a interface.",err,error,*999)
            ENDIF
            IF(data_points_region_dimensions==mesh_region_dimensions) THEN !Dimension has to be equal
              ALLOCATE(data_projection,stat=err)
              IF(err/=0) CALL flag_error("Could not allocate data projection.",err,error,*999)
              data_projection%USER_NUMBER=data_projection_user_number
              data_projection%LABEL=""
              CALL tree_item_insert(data_points%DATA_PROJECTIONS_TREE,data_projection_user_number, &
                & data_points%NUMBER_OF_DATA_PROJECTIONS+1,insert_status,err,error,*999)
              data_projection%DATA_PROJECTION_FINISHED=.false.
              data_projection%DATA_POINTS=>data_points
              data_projection%MESH=>mesh
              data_projection%COORDINATE_SYSTEM_DIMENSIONS=data_points_region_dimensions
              data_projection%MAXIMUM_ITERATION_UPDATE=0.5_dp
              data_projection%MAXIMUM_NUMBER_OF_ITERATIONS=25
              !Default always project to boundaries faces/lines when mesh dimension is equal to region dimension. If mesh dimension is less, project to all elements
              IF(mesh%NUMBER_OF_DIMENSIONS<data_points_region_dimensions) THEN !mesh dimension < data dimension
                data_projection%NUMBER_OF_XI=mesh%NUMBER_OF_DIMENSIONS
                data_projection%PROJECTION_TYPE=data_projection_all_elements_projection_type
              ELSE
                SELECT CASE(mesh%NUMBER_OF_DIMENSIONS) !mesh dimension = data dimension
                CASE (2)
                  data_projection%NUMBER_OF_XI=1
                  data_projection%PROJECTION_TYPE=data_projection_boundary_lines_projection_type
                CASE (3)
                  data_projection%NUMBER_OF_XI=2
                  data_projection%PROJECTION_TYPE=data_projection_boundary_faces_projection_type
                CASE DEFAULT
                  CALL flag_error("Mesh dimensions out of bond [1,3].",err,error,*999)
                END SELECT
              ENDIF
              SELECT CASE(data_projection%NUMBER_OF_XI) !mesh dimension = data dimension
              CASE (1)
                data_projection%NUMBER_OF_CLOSEST_ELEMENTS=2
              CASE (2)
                data_projection%NUMBER_OF_CLOSEST_ELEMENTS=4
              CASE (3)
                data_projection%NUMBER_OF_CLOSEST_ELEMENTS=8
              CASE DEFAULT
                CALL flag_error("Mesh dimensions out of bond [1,3].",err,error,*999)
              END SELECT
              ALLOCATE(data_projection%STARTING_XI(data_projection%NUMBER_OF_XI),stat=err)
              IF(err/=0) CALL flag_error("Could not allocate data points data projection starting xi.",err,error,*999)
              DO xi_idx=1,data_projection%NUMBER_OF_XI
                data_projection%STARTING_XI(xi_idx)=0.5_dp 
              ENDDO !xi_idx
              data_projection%ABSOLUTE_TOLERANCE=1.0e-8_dp
              data_projection%RELATIVE_TOLERANCE=1.0e-6_dp
              IF(data_points%NUMBER_OF_DATA_PROJECTIONS>0) THEN
                ALLOCATE(new_data_projections_ptr(data_points%NUMBER_OF_DATA_PROJECTIONS+1),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate new data projections.",err,error,*999)
                DO data_projection_idx=1,data_points%NUMBER_OF_DATA_PROJECTIONS
                  new_data_projections_ptr(data_projection_idx)%PTR=>data_points%DATA_PROJECTIONS(data_projection_idx)%PTR
                ENDDO !xi_idx
              ELSE IF(data_points%NUMBER_OF_DATA_PROJECTIONS==0) THEN
                ALLOCATE(new_data_projections_ptr(data_points%NUMBER_OF_DATA_PROJECTIONS+1),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate new data projections.",err,error,*999)
              ELSE
                CALL flag_error("The number of data projections is < 0.",err,error,*999)
              ENDIF
              !Return the pointer
              new_data_projections_ptr(data_points%NUMBER_OF_DATA_PROJECTIONS+1)%PTR=>data_projection
              CALL move_alloc(new_data_projections_ptr,data_points%DATA_PROJECTIONS)
              data_points%NUMBER_OF_DATA_PROJECTIONS=data_points%NUMBER_OF_DATA_PROJECTIONS+1
              data_projection%GLOBAL_NUMBER=data_points%NUMBER_OF_DATA_PROJECTIONS
            ELSE
              CALL flag_error("Dimensions bewtween the mesh region/interface and data points region/interface does not match.", &
                & err,error,*999)
            ENDIF !DATA_POINTS_REGION_DIMENSIONS==MESH_REGION_DIMENSIONS
          ENDIF !ASSOCIATED(DATA_PROJECTION)
        ELSE
          CALL flag_error("Mesh is not associated.",err,error,*999)
        ENDIF !ASSOCIATED(MESH)
      ELSE
        CALL flag_error("Data points have not been finished.",err,error,*999)
      ENDIF !DATA_POINTS%DATA_POINTS_FINISHED
    ELSE
      CALL flag_error("Data points is not associated.",err,error,*999)
    ENDIF !ASSOCIATED(DATA_POINTS)

    exits("DATA_PROJECTION_CREATE_START_DATA_POINTS")
    RETURN
999 errorsexits("DATA_PROJECTION_CREATE_START_DATA_POINTS",err,error)
    RETURN 1

  END SUBROUTINE data_projection_create_start_data_points

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

  SUBROUTINE dataprojection_datapointcheckexist(DataProjection,dataPointUserNumber,dataPointExist,dataPointGlobalNumber,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DataProjection
    INTEGER(INTG) :: dataPointUserNumber
    LOGICAL, INTENT(OUT) :: dataPointExist
    INTEGER(INTG), INTENT(OUT) :: dataPointGlobalNumber
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    TYPE(data_points_type), POINTER :: dataPoints
    TYPE(tree_node_type), POINTER :: treeNode

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

    datapointexist=.false.
    datapointglobalnumber=0
    IF(ASSOCIATED(dataprojection)) THEN
      datapoints=>dataprojection%DATA_POINTS
      IF(ASSOCIATED(datapoints)) THEN
        NULLIFY(treenode)
        CALL tree_search(datapoints%DATA_POINTS_TREE,datapointusernumber,treenode,err,error,*999)
        IF(ASSOCIATED(treenode)) THEN
          CALL tree_node_value_get(datapoints%DATA_POINTS_TREE,treenode,datapointglobalnumber,err,error,*999)
          datapointexist=.true.
        ENDIF
      ELSE
        CALL flag_error("Data points is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_DataPointCheckExist")
    RETURN
999 errorsexits("DataProjection_DataPointCheckExist",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_datapointcheckexist

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

  SUBROUTINE dataprojection_dataprojectionresultinitialise(DATA_PROJECTION,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: NUMBER_OF_DATA_POINTS
    INTEGER(INTG) :: data_point_idx

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

    !NULLIFY(PROJECTED_POINTS)
    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(ASSOCIATED(data_projection%DATA_POINTS)) THEN
          IF(data_projection%DATA_POINTS%DATA_POINTS_FINISHED) THEN
            number_of_data_points = data_projection%DATA_POINTS%NUMBER_OF_DATA_POINTS
            ALLOCATE(data_projection%DATA_PROJECTION_RESULTS(number_of_data_points),stat=err)
            IF(err/=0) CALL flag_error("Could not allocate data projection data projection results.",err,error,*999)
            DO data_point_idx=1,number_of_data_points
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%USER_NUMBER=data_projection%DATA_POINTS%DATA_POINTS( &
                & data_point_idx)%USER_NUMBER
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%DISTANCE=0.0_dp
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%ELEMENT_NUMBER=0
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%ELEMENT_FACE_NUMBER=0
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%ELEMENT_LINE_NUMBER=0
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%EXIT_TAG=data_projection_exit_tag_no_element
              ALLOCATE(data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI(data_projection%NUMBER_OF_XI),stat=err)
        ALLOCATE(data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector( &
                & data_projection%COORDINATE_SYSTEM_DIMENSIONS),stat=err)
              IF(err/=0) CALL flag_error("Could not allocate data projection data projection results "// &
                & "("//trim(number_to_vstring(data_point_idx,"*",err,error))//") xi.",err,error,*999)
              data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI(1:data_projection%NUMBER_OF_XI)= &
                & data_projection%STARTING_XI(1:data_projection%NUMBER_OF_XI)
        data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector( &
                & 1:data_projection%COORDINATE_SYSTEM_DIMENSIONS)=0.0_dp
            ENDDO !data_point_idx
          ELSE
            CALL flag_error("Data projection data points have not been finished.",err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection data points is not associated.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_DataProjectionResultInitialise")
    RETURN
999 errors("DataProjection_DataProjectionResultInitialise",err,error)
    exits("DataProjection_DataProjectionResultInitialise")
    RETURN 1

  END SUBROUTINE dataprojection_dataprojectionresultinitialise

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

  SUBROUTINE data_projection_destroy(DATA_PROJECTION,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      CALL data_projection_finalise(data_projection,err,error,*999)
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_DESTROY")
    RETURN
999 errorsexits("DATA_PROJECTION_DESTROY",err,error)
    RETURN 1

  END SUBROUTINE data_projection_destroy

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

  SUBROUTINE data_projection_finalise(DATA_PROJECTION,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: dataPointIdx

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(ALLOCATED(data_projection%STARTING_XI)) THEN
        DEALLOCATE(data_projection%STARTING_XI)
      ENDIF
      IF(ALLOCATED(data_projection%DATA_PROJECTION_RESULTS)) THEN
        DO datapointidx=1,SIZE(data_projection%DATA_PROJECTION_RESULTS,1)
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%USER_NUMBER=0
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%DISTANCE=0
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%ELEMENT_NUMBER=0
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%ELEMENT_FACE_NUMBER=0
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%ELEMENT_LINE_NUMBER=0
          data_projection%DATA_PROJECTION_RESULTS(datapointidx)%EXIT_TAG=0
          IF(ALLOCATED(data_projection%DATA_PROJECTION_RESULTS(datapointidx)%XI)) THEN
            DEALLOCATE(data_projection%DATA_PROJECTION_RESULTS(datapointidx)%XI)
          ENDIF
          IF(ALLOCATED(data_projection%DATA_PROJECTION_RESULTS(datapointidx)%projectionVector)) THEN
            DEALLOCATE(data_projection%DATA_PROJECTION_RESULTS(datapointidx)%projectionVector)
          ENDIF
        ENDDO !dataPointIdx
        DEALLOCATE(data_projection%DATA_PROJECTION_RESULTS)
      ENDIF
      DEALLOCATE(data_projection)
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_FINALISE")
    RETURN
999 errorsexits("DATA_PROJECTION_FINALISE",err,error)
    RETURN 1

  END SUBROUTINE data_projection_finalise

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

  SUBROUTINE dataprojection_datapointglobalnumberget(DATA_PROJECTION,USER_NUMBER,GLOBAL_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: GLOBAL_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(data_points_type), POINTER  :: DATA_POINTS
    TYPE(varying_string) :: LOCAL_ERROR
    TYPE(tree_node_type), POINTER :: TREE_NODE

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

    IF(ASSOCIATED(data_projection)) THEN
      data_points=>data_projection%DATA_POINTS
      IF(ASSOCIATED(data_points)) THEN
        NULLIFY(tree_node)
        CALL tree_search(data_points%DATA_POINTS_TREE,user_number,tree_node,err,error,*999)
        IF(ASSOCIATED(tree_node)) THEN
          CALL tree_node_value_get(data_points%DATA_POINTS_TREE,tree_node,global_number,err,error,*999)
        ELSE
          local_error="Tree node is not associates (cannot find the user number "//trim(number_to_vstring(user_number,"*",err, &
            & error))//"."
          CALL flag_error(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data points is not associated.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_DataPointGlobalNumberGet")
    RETURN
999 errorsexits("DataProjection_DataPointGlobalNumberGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_datapointglobalnumberget

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


  SUBROUTINE dataprojection_datapointsprojectionevaluate(DATA_PROJECTION,PROJECTION_FIELD,ERR,ERROR,*) !optimising

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_type), POINTER :: PROJECTION_FIELD
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    TYPE(field_interpolation_parameters_ptr_type), POINTER :: INTERPOLATION_PARAMETERS(:)
    TYPE(field_interpolated_point_ptr_type), POINTER :: INTERPOLATED_POINTS(:)
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    TYPE(data_points_type), POINTER :: DATA_POINTS
    TYPE(decomposition_type), POINTER :: DECOMPOSITION
    TYPE(domain_type), POINTER :: DOMAIN
    !TYPE(DATA_PROJECTION_PROJECTED_POINTS_TYPE), POINTER :: PROJECTED_POINTS
    INTEGER(INTG) :: MY_COMPUTATIONAL_NODE,NUMBER_COMPUTATIONAL_NODES
    INTEGER(INTG) :: MESH_COMPONENT_NUMBER,NUMBER_OF_DATA_POINTS
    INTEGER(INTG) :: NUMBER_OF_ELEMENTS,NUMBER_OF_FACES,NUMBER_OF_LINES,NUMBER_OF_CANDIDATES
    INTEGER(INTG) :: NUMBER_OF_CLOSEST_CANDIDATES,TOTAL_NUMBER_OF_CLOSEST_CANDIDATES,REDUCED_NUMBER_OF_CLOSEST_CANDIDATES
    INTEGER(INTG), ALLOCATABLE :: GLOBAL_TO_LOCAL_NUMBER_OF_CLOSEST_CANDIDATES(:)
    INTEGER(INTG), ALLOCATABLE :: CANDIDATE_ELEMENTS(:),CLOSEST_ELEMENTS(:,:),CANDIDATE_FACES(:),CLOSEST_FACES(:,:)
    !INTEGER(INTG) :: NUMBER_OF_CANDIDATE_LINES
    REAL(DP), ALLOCATABLE :: CLOSEST_DISTANCES(:,:),GLOBAL_CLOSEST_DISTANCES(:,:)
    INTEGER(INTG), ALLOCATABLE :: GLOBAL_NUMBER_OF_CLOSEST_CANDIDATES(:)
    INTEGER(INTG), ALLOCATABLE :: GLOBAL_MPI_DISPLACEMENTS(:),SORTING_IND_1(:),SORTING_IND_2(:)
    INTEGER(INTG), ALLOCATABLE :: GLOBAL_NUMBER_OF_PROJECTED_POINTS(:)
    INTEGER(INTG) :: MPI_CLOSEST_DISTANCES,DATA_PROJECTION_GLOBAL_NUMBER
    INTEGER(INTG) :: MPI_IERROR
    !LOGICAL :: ELEMENT_FOUND
    !LOGICAL, ALLOCATABLE :: PROJECTION_CONVERGED(:)
    INTEGER(INTG), ALLOCATABLE :: PROJECTED_ELEMENT(:),PROJECTED_FACE(:),PROJECTION_EXIT_TAG(:)
    REAL(DP), ALLOCATABLE :: PROJECTED_DISTANCE(:,:),PROJECTED_XI(:,:), PROJECTION_VECTORS(:,:)

    INTEGER(INTG) :: ne,nse,ncn,ni,localElementNumber
    INTEGER(INTG) :: temp_number,start_idx,finish_idx,data_point_idx

    LOGICAL :: BOUNDARY_PROJECTION,elementExists,ghostElement

    !INTEGER(INTG) :: NUMBER_OF_PROJECTED_POINTS,NUMBER_OF_PROJECTED_POINTS_2

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

    NULLIFY(interpolation_parameters)
    NULLIFY(interpolated_points)
    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(ASSOCIATED(projection_field)) THEN
          IF(projection_field%FIELD_FINISHED) THEN
            IF(ASSOCIATED(data_projection%MESH,projection_field%DECOMPOSITION%MESH)) THEN
              data_projection%PROJECTION_FIELD=>projection_field
              data_projection_global_number=data_projection%GLOBAL_NUMBER
              data_points=>data_projection%DATA_POINTS
              CALL field_interpolation_parameters_initialise(projection_field,interpolation_parameters,err,error,*998, &
                & field_geometric_components_type)
              CALL field_interpolated_points_initialise(interpolation_parameters,interpolated_points,err,error,*998, &
                & field_geometric_components_type)
              interpolated_point=>interpolated_points(field_u_variable_type)%PTR
              decomposition=>projection_field%DECOMPOSITION
              mesh_component_number=projection_field%DECOMPOSITION%MESH_COMPONENT_NUMBER
              domain=>projection_field%DECOMPOSITION%DOMAIN(mesh_component_number)%PTR
              number_of_data_points=data_points%NUMBER_OF_DATA_POINTS
              my_computational_node=computational_node_number_get(err,error)
              number_computational_nodes=computational_environment%NUMBER_COMPUTATIONAL_NODES
              boundary_projection=(data_projection%PROJECTION_TYPE==data_projection_boundary_lines_projection_type).OR.( &
                & data_projection%PROJECTION_TYPE==data_projection_boundary_faces_projection_type)
              !#####################################################################################################################
              !find elements/faces/lines inside current computational node, get the boundary faces/lines only if asked
              !the elements/faces/lines are required to perform projection of points in the current computational node
              !the are all pre-allocated to the maximum array length (i.e. NUMBER_OF_ELEMENTS), but only up to the NUMBER_OF_CANDIDATESth index are assigned
              number_of_elements=domain%TOPOLOGY%ELEMENTS%NUMBER_OF_ELEMENTS
              number_of_faces=domain%TOPOLOGY%FACES%NUMBER_OF_FACES
              number_of_lines=domain%TOPOLOGY%LINES%NUMBER_OF_LINES
              number_of_candidates=0
              IF(ALLOCATED(data_projection%candidateElementNumbers) .AND. ALLOCATED(data_projection%localFaceLineNumbers)) THEN
                SELECT CASE(data_projection%PROJECTION_TYPE)
                CASE (data_projection_boundary_lines_projection_type)  !identify all non-ghost boundary lines
                  ALLOCATE(candidate_elements(number_of_lines),stat=err)
                  IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                  ALLOCATE(candidate_faces(number_of_lines),stat=err)
                  IF(err/=0) CALL flag_error("Could not allocate candidate lines.",err,error,*999)
                  DO ne=1,SIZE(data_projection%candidateElementNumbers,1) !Loop through all candidate element defined by user number
                    CALL decomposition_topology_element_check_exists(decomposition%TOPOLOGY,data_projection% &
                      & candidateelementnumbers(ne),elementexists,localelementnumber,ghostelement,err,error,*999) !Check if element exists on current domain, get local number
                    IF((elementexists) .AND. (.NOT.ghostelement)) THEN !Get non-ghost elements
                      number_of_candidates=number_of_candidates+1
                      candidate_elements(number_of_candidates)=localelementnumber
                      candidate_faces(number_of_candidates)=data_projection%localFaceLineNumbers(ne) !Store element line number for line projection type
                    ENDIF
                  ENDDO !ne
                CASE (data_projection_boundary_faces_projection_type) !identify all non-ghost boundary faces
                  IF(decomposition%MESH%NUMBER_OF_DIMENSIONS>=2) THEN
                    ALLOCATE(candidate_elements(number_of_faces),stat=err)
                    IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                    ALLOCATE(candidate_faces(number_of_faces),stat=err)
                    IF(err/=0) CALL flag_error("Could not allocate candidate faces.",err,error,*999)
                    DO ne=1,SIZE(data_projection%candidateElementNumbers,1) !Loop through all candidate element defined by user number
                      CALL decomposition_topology_element_check_exists(decomposition%TOPOLOGY,data_projection% &
                        & candidateelementnumbers(ne),elementexists,localelementnumber,ghostelement,err,error,*999) !Check if element exists on current domain, get local number
                      IF((elementexists) .AND. (.NOT.ghostelement)) THEN !Get non-ghost elements
                        number_of_candidates=number_of_candidates+1
                        candidate_elements(number_of_candidates)=localelementnumber
                        candidate_faces(number_of_candidates)=data_projection%localFaceLineNumbers(ne) !Store element face number for face projection type
                      ENDIF
                    ENDDO !ne
                  ELSE
                    CALL flag_error("Decomposition mesh number of dimensions has to be 2 or greater.",err,error,*999)
                  ENDIF
                CASE (data_projection_all_elements_projection_type) !identify all non-ghost elements
                  IF(data_projection%NUMBER_OF_XI==decomposition%MESH%NUMBER_OF_DIMENSIONS) THEN
                    ALLOCATE(candidate_elements(number_of_elements),stat=err)
                    IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                    DO ne=1,SIZE(data_projection%candidateElementNumbers,1) !Loop through all candidate element defined by user number
                      CALL decomposition_topology_element_check_exists(decomposition%TOPOLOGY,data_projection% &
                        & candidateelementnumbers(ne),elementexists,localelementnumber,ghostelement,err,error,*999) !Check if element exists on current domain, get local number
                      IF((elementexists) .AND. (.NOT.ghostelement)) THEN !Get non-ghost elements
                        number_of_candidates=number_of_candidates+1
                        candidate_elements(number_of_candidates)=localelementnumber
                      ENDIF
                    ENDDO !ne
                  ELSE
                    CALL flag_error("Data projection number of xi has to equal to decomposition mesh number of dimensions",err, &
                      & error,*999)
                  ENDIF
                CASE DEFAULT
                  CALL flag_error("No match for data projection type found",err,error,*999)
                END SELECT !DATA_PROJECTION%PROJECTION_TYPE
              ELSE !If user didn't define candidate element number
                SELECT CASE(data_projection%PROJECTION_TYPE)
                  CASE (data_projection_boundary_lines_projection_type)  !identify all non-ghost boundary lines
                    ALLOCATE(candidate_elements(number_of_lines),stat=err)
                    IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                    ALLOCATE(candidate_faces(number_of_lines),stat=err)
                    IF(err/=0) CALL flag_error("Could not allocate candidate lines.",err,error,*999)
                    DO ne=1,domain%MAPPINGS%ELEMENTS%NUMBER_OF_LOCAL
                      IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BOUNDARY_ELEMENT) THEN
                        DO nse=1,SIZE(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%ELEMENT_LINES,1)
                          temp_number=decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%ELEMENT_LINES(nse)
                          IF(decomposition%TOPOLOGY%LINES%LINES(temp_number)%BOUNDARY_LINE) THEN
                            number_of_candidates=number_of_candidates+1
                            candidate_faces(number_of_candidates)=nse
                            candidate_elements(number_of_candidates)=ne
                          ENDIF !boundary line
                        ENDDO !nse
                      ENDIF !boundary element
                    ENDDO !ne
                  CASE (data_projection_boundary_faces_projection_type) !identify all non-ghost boundary faces
                    IF(decomposition%MESH%NUMBER_OF_DIMENSIONS>=2) THEN
                      ALLOCATE(candidate_elements(number_of_faces),stat=err)
                      IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                      ALLOCATE(candidate_faces(number_of_faces),stat=err)
                      IF(err/=0) CALL flag_error("Could not allocate candidate faces.",err,error,*999)
                      DO ne=1,domain%MAPPINGS%ELEMENTS%NUMBER_OF_LOCAL
                        IF(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%BOUNDARY_ELEMENT) THEN
                          DO nse=1,SIZE(decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%ELEMENT_FACES,1)
                            temp_number=decomposition%TOPOLOGY%ELEMENTS%ELEMENTS(ne)%ELEMENT_FACES(nse)
                            IF(decomposition%TOPOLOGY%FACES%FACES(temp_number)%BOUNDARY_FACE) THEN
                              number_of_candidates=number_of_candidates+1
                              candidate_faces(number_of_candidates)=nse
                              candidate_elements(number_of_candidates)=ne
                            ENDIF !boundary face
                          ENDDO !nse
                        ENDIF !boundary element
                      ENDDO !ne
                    ELSE
                      CALL flag_error("Decomposition mesh number of dimensions has to be 2 or greater.",err,error,*999)
                    ENDIF
                  CASE (data_projection_all_elements_projection_type) !identify all non-ghost elements
                    IF(data_projection%NUMBER_OF_XI==decomposition%MESH%NUMBER_OF_DIMENSIONS) THEN
                      ALLOCATE(candidate_elements(number_of_elements),stat=err)
                      IF(err/=0) CALL flag_error("Could not allocate candidate elements.",err,error,*999)
                      DO ne=1,domain%MAPPINGS%ELEMENTS%NUMBER_OF_LOCAL
                        number_of_candidates=number_of_candidates+1
                        candidate_elements(number_of_candidates)=ne
                      ENDDO
                    ELSE
                      CALL flag_error("Data projection number of xi has to equal to decomposition mesh number of dimensions",err, &
                        & error,*999)
                    ENDIF
                  CASE DEFAULT
                    CALL flag_error("No match for data projection type found",err,error,*999)
                END SELECT !DATA_PROJECTION%PROJECTION_TYPE
              ENDIF
              !#####################################################################################################################
              !find the clostest elements/faces/lines for each point in the current computational node base on starting xi
              !the clostest elements/faces/lines are required to shrink down on the list of possible projection candiates
              number_of_closest_candidates=min(data_projection%NUMBER_OF_CLOSEST_ELEMENTS,number_of_candidates)
              ALLOCATE(closest_elements(number_of_data_points,number_of_closest_candidates),stat=err)!the information for each data point has to be stored in the corresponding rows for them to be contiguous in memory for easy MPI access
              IF(err/=0) CALL flag_error("Could not allocate closest elements.",err,error,*999)
              IF(boundary_projection) THEN
                ALLOCATE(closest_faces(number_of_data_points,number_of_closest_candidates),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate closest sub element.",err,error,*999)
              ENDIF
              ALLOCATE(closest_distances(number_of_data_points,number_of_closest_candidates),stat=err)!the information for each data point has to be stored in the corresponding rows for them to be contiguous in memory for easy MPI access
              IF(err/=0) CALL flag_error("Could not allocate closest distances.",err,error,*999)
              SELECT CASE(data_projection%PROJECTION_TYPE)
                CASE (data_projection_boundary_lines_projection_type) !find closest candidate lines
                  DO data_point_idx=1,number_of_data_points
                    CALL data_projection_closest_lines_find(data_projection,interpolated_point, &
                      & data_points%DATA_POINTS(data_point_idx)%position,candidate_elements, &
                      & candidate_faces,number_of_candidates,closest_elements(data_point_idx,:),closest_faces( &
                      & data_point_idx,:),closest_distances(data_point_idx,:),err,error,*999)
                  ENDDO !data_point_idx
                CASE (data_projection_boundary_faces_projection_type) !find closest candidate faces
                  DO data_point_idx=1,number_of_data_points
                    CALL data_projection_closest_faces_find(data_projection,interpolated_point, &
                      & data_points%DATA_POINTS(data_point_idx)%position,candidate_elements, &
                      & candidate_faces,number_of_candidates,closest_elements(data_point_idx,:),closest_faces( &
                      & data_point_idx,:),closest_distances(data_point_idx,:),err,error,*999)
                  ENDDO !data_point_idx
                CASE (data_projection_all_elements_projection_type) !find closest candidate elements
                  DO data_point_idx=1,number_of_data_points
                    CALL data_projection_closest_elements_find(data_projection,interpolated_point,data_points%DATA_POINTS( &
                      & data_point_idx)%position,candidate_elements,number_of_candidates,closest_elements(data_point_idx,:), &
                      & closest_distances(data_point_idx,:),err,error,*999)
                    !CLOSEST_ELEMENTS(data_point_idx,:)=DOMAIN%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(CLOSEST_ELEMENTS(data_point_idx,:)) !local to global element number mapping
                  ENDDO !data_point_idx
                CASE DEFAULT
                  CALL flag_error("No match for data projection type found",err,error,*999)
              END SELECT
              !#####################################################################################################################
              !Newton project data point to the list of closest elements, faces or lines
              !project the data points to each of the closest elements, use MPI if number of computational nodes is greater than 1
              IF(number_computational_nodes>1) THEN !use mpi
                !allocate arrays for mpi communication
                ALLOCATE(global_to_local_number_of_closest_candidates(number_of_data_points),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate global to local number of closest elements.",err,error,*999)
                ALLOCATE(global_number_of_closest_candidates(number_computational_nodes),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate all number of closest candidates.",err,error,*999)
                ALLOCATE(global_mpi_displacements(number_computational_nodes),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate all displacements.",err,error,*999)
                ALLOCATE(global_number_of_projected_points(number_computational_nodes),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate all number of projected points.",err,error,*999)
                ALLOCATE(projection_exit_tag(number_of_data_points),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate projected.",err,error,*999)
                ALLOCATE(projected_element(number_of_data_points),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate projected element.",err,error,*999)
                IF(boundary_projection) THEN
                  ALLOCATE(projected_face(number_of_data_points),stat=err)
                  IF(err/=0) CALL flag_error("Could not allocate projected sub element.",err,error,*999)
                ENDIF
                ALLOCATE(projected_distance(2,number_of_data_points),stat=err) !PROJECTED_DISTANCE(2,:) stores the compuational node number, the information for each data point has to be stored in the corresponding column for MPI_ALLREDUCE with location return to work
                IF(err/=0) CALL flag_error("Could not allocate projected distance.",err,error,*999)
                ALLOCATE(projected_xi(data_projection%NUMBER_OF_XI,number_of_data_points),stat=err) !the information for each data point is stored in the corresponding column to be consistent with PROJECTED_DISTANCE
                ALLOCATE(projection_vectors(data_projection%COORDINATE_SYSTEM_DIMENSIONS, number_of_data_points), stat=err) !ZJW: The projection vector for each data point is stored in the corresponding colum to be consistent with PROJECTED_DISTANCE.
                IF(err/=0) CALL flag_error("Could not allocate projected.",err,error,*999)
                ALLOCATE(sorting_ind_2(number_of_data_points),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate sorting ind 2.",err,error,*999)
                !gather and distribute the number of closest elements from all computational nodes
                CALL mpi_allgather(number_of_closest_candidates,1,mpi_integer,global_number_of_closest_candidates,1,mpi_integer, &
                  & computational_environment%MPI_COMM,mpi_ierror)
                CALL mpi_error_check("MPI_ALLGATHER",mpi_ierror,err,error,*999)
                total_number_of_closest_candidates=sum(global_number_of_closest_candidates,1) !sum of all number of closest candidates from all computational nodes
                !allocate arrays to store information gathered from all computational node
                ALLOCATE(global_closest_distances(number_of_data_points,total_number_of_closest_candidates),stat=err) !the information for each data point is stored in the corresponding row so they are contiguous in memory for easy MPI access
                IF(err/=0) CALL flag_error("Could not allocate all closest distances.",err,error,*999)
                ALLOCATE(sorting_ind_1(total_number_of_closest_candidates),stat=err)
                IF(err/=0) CALL flag_error("Could not allocate sorting ind 1.",err,error,*999)
                !MPI:create and commit MPI_TYPE_CONTIGUOUS
                CALL mpi_type_contiguous(number_of_data_points,mpi_double_precision,mpi_closest_distances,mpi_ierror)
                CALL mpi_error_check("MPI_TYPE_CONTIGUOUS",mpi_ierror,err,error,*999)
                CALL mpi_type_commit(mpi_closest_distances,mpi_ierror)
                CALL mpi_error_check("MPI_TYPE_COMMIT",mpi_ierror,err,error,*999)
                !create displacement vectors for MPI_ALLGATHERV
                global_mpi_displacements(1)=0
                DO ncn=1,(number_computational_nodes-1)
                  global_mpi_displacements(ncn+1)=global_mpi_displacements(ncn)+global_number_of_closest_candidates(ncn)
                ENDDO
                !MPI:shares closest element distances between all domains
                CALL mpi_allgatherv(closest_distances(1,1),number_of_closest_candidates,mpi_closest_distances, &
                  & global_closest_distances,global_number_of_closest_candidates,global_mpi_displacements,mpi_closest_distances, &
                  & computational_environment%MPI_COMM,mpi_ierror)
                CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
                reduced_number_of_closest_candidates=min(data_projection%NUMBER_OF_CLOSEST_ELEMENTS, &
                  & total_number_of_closest_candidates)
                projected_distance(2,:)=my_computational_node
                ! find the globally closest distances in the current domain
                DO data_point_idx=1,number_of_data_points
                  CALL bubble_isort(global_closest_distances(data_point_idx,:),sorting_ind_1,err,error,*999)
                  sorting_ind_1(1:total_number_of_closest_candidates)=sorting_ind_1(1:total_number_of_closest_candidates)- &
                    & global_mpi_displacements(my_computational_node+1) !shift the index to current computational node
                  global_to_local_number_of_closest_candidates(data_point_idx)=0
                  DO ne=1,reduced_number_of_closest_candidates
                    !sorted index indicates it is in the current computational domain
                    IF((sorting_ind_1(ne)>=1).AND.(sorting_ind_1(ne)<=global_number_of_closest_candidates(my_computational_node &
                      & +1)))global_to_local_number_of_closest_candidates(data_point_idx)= &
                      & global_to_local_number_of_closest_candidates(data_point_idx)+1
                  ENDDO
                  projected_distance(1,data_point_idx)=global_closest_distances(data_point_idx,total_number_of_closest_candidates) !assign initial distance to something large
                ENDDO
                SELECT CASE(data_projection%PROJECTION_TYPE)
                  CASE (data_projection_boundary_lines_projection_type) !Newton project to closest lines, and find miminum projection
                    DO data_point_idx=1,number_of_data_points
                      number_of_closest_candidates=global_to_local_number_of_closest_candidates(data_point_idx)
                      IF(number_of_closest_candidates>0) THEN
                        CALL data_projection_newton_lines_evaluate(data_projection,interpolated_point, &
                          & data_points%DATA_POINTS(data_point_idx)%position,closest_elements( &
                          & data_point_idx,1:number_of_closest_candidates),closest_faces(data_point_idx,1: &
                          & number_of_closest_candidates),projection_exit_tag(data_point_idx),projected_element(data_point_idx),  &
                          & projected_face(data_point_idx),projected_distance(1,data_point_idx),projected_xi(:,data_point_idx), &
                          & projection_vectors(:, data_point_idx),err,error,*999)
                        projected_element(data_point_idx)=domain%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(projected_element( &
                          & data_point_idx)) !map the element number to global number
                      ENDIF
                    ENDDO
                  CASE (data_projection_boundary_faces_projection_type) !find closest candidate faces
                    DO data_point_idx=1,number_of_data_points
                      number_of_closest_candidates=global_to_local_number_of_closest_candidates(data_point_idx)
                      IF(number_of_closest_candidates>0) THEN
                        CALL data_projection_newton_faces_evaluate(data_projection,interpolated_point, &
                          & data_points%DATA_POINTS(data_point_idx)%position,closest_elements( &
                          & data_point_idx,1:number_of_closest_candidates),closest_faces(data_point_idx, &
                          & 1:number_of_closest_candidates),projection_exit_tag(data_point_idx),projected_element(data_point_idx), &
                          & projected_face(data_point_idx),projected_distance(1,data_point_idx),projected_xi(:,data_point_idx), &
                          & projection_vectors(:, data_point_idx),err,error,*999)
                        projected_element(data_point_idx)=domain%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(projected_element( &
                          & data_point_idx)) !map the element number to global number
                      ENDIF
                    ENDDO
                  CASE (data_projection_all_elements_projection_type) !find closest candidate elements
                    SELECT CASE(data_projection%NUMBER_OF_XI)
                      CASE (1) !1D element
                        DO data_point_idx=1,number_of_data_points
                          number_of_closest_candidates=global_to_local_number_of_closest_candidates(data_point_idx)
                          IF(number_of_closest_candidates>0) THEN
                            CALL data_projection_newton_elements_evaluate_1(data_projection,interpolated_point,data_points% &
                              & data_points(data_point_idx)%position,closest_elements(data_point_idx, &
                              & 1:number_of_closest_candidates),projection_exit_tag(data_point_idx), &
                              & projected_element(data_point_idx),projected_distance(1,data_point_idx), &
                              & projected_xi(:,data_point_idx),projection_vectors(:, data_point_idx),err,error,*999)
                            projected_element(data_point_idx)=domain%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(projected_element( &
                              & data_point_idx)) !map the element number to global number

                          ENDIF
                        ENDDO
                      CASE (2) !2D element
                        DO data_point_idx=1,number_of_data_points
                          number_of_closest_candidates=global_to_local_number_of_closest_candidates(data_point_idx)
                          IF(number_of_closest_candidates>0) THEN
                            CALL data_projection_newton_elements_evaluate_2(data_projection,interpolated_point,data_points% &
                              & data_points(data_point_idx)%position,closest_elements(data_point_idx, &
                              & 1:number_of_closest_candidates),projection_exit_tag(data_point_idx), &
                              & projected_element(data_point_idx),projected_distance(1,data_point_idx), &
                              & projected_xi(:,data_point_idx),projection_vectors(:, data_point_idx), &
            & err,error,*999)
                            projected_element(data_point_idx)=domain%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(projected_element( &
                              & data_point_idx)) !map the element number to global number
                          ENDIF
                        ENDDO
                      CASE (3) !3D element
                        DO data_point_idx=1,number_of_data_points
                          number_of_closest_candidates=global_to_local_number_of_closest_candidates(data_point_idx)
                          IF(number_of_closest_candidates>0) THEN
                            CALL data_projection_newton_elements_evaluate_3(data_projection,interpolated_point,data_points% &
                              & data_points(data_point_idx)%position,closest_elements(data_point_idx, &
                              & 1:number_of_closest_candidates),projection_exit_tag(data_point_idx), &
                              & projected_element(data_point_idx),projected_distance(1,data_point_idx), &
                              & projected_xi(:,data_point_idx),projection_vectors(:, data_point_idx),err,error,*999)
                            projected_element(data_point_idx)=domain%MAPPINGS%ELEMENTS%LOCAL_TO_GLOBAL_MAP(projected_element( &
                              & data_point_idx)) !map the element number to global number
                          ENDIF
                        ENDDO
                      CASE DEFAULT
                        CALL flag_error("Data projection number of xi is invalid",err,error,*999)
                    END SELECT
                  CASE DEFAULT
                    CALL flag_error("No match for data projection type found",err,error,*999)
                END SELECT
                !MPI:find the shortest projected distance in all domains
                CALL mpi_allreduce(mpi_in_place,projected_distance,number_of_data_points,mpi_2double_precision,mpi_minloc, &
                  & computational_environment%MPI_COMM,mpi_ierror)
                CALL mpi_error_check("MPI_ALLREDUCE",mpi_ierror,err,error,*999)
                !sort the computational node/rank from 0 to number of computational node
                CALL bubble_isort(projected_distance(2,:),sorting_ind_2,err,error,*999)
                DO ncn=0,(number_computational_nodes-1)
                  global_number_of_projected_points(ncn+1)=count(abs(projected_distance(2,:)-REAL(ncn))<zero_tolerance)
                ENDDO !ncn
                start_idx=sum(global_number_of_projected_points(1:my_computational_node))+1
                finish_idx=start_idx+global_number_of_projected_points(my_computational_node+1)-1
                !create displacement vectors for MPI_ALLGATHERV
                DO ncn=1,(number_computational_nodes-1)
                  global_mpi_displacements(ncn+1)=global_mpi_displacements(ncn)+global_number_of_projected_points(ncn)
                ENDDO !ncn
                !MPI:shares minimum projection informationg bewteen all domains
                CALL mpi_allgatherv(projected_element(sorting_ind_2(start_idx:finish_idx)),global_number_of_projected_points( &
                  & my_computational_node+1),mpi_integer,projected_element,global_number_of_projected_points, &
                  & global_mpi_displacements,mpi_integer,computational_environment%MPI_COMM,mpi_ierror) !PROJECTED_ELEMENT
                CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
                IF(boundary_projection) THEN
                  CALL mpi_allgatherv(projected_face(sorting_ind_2(start_idx:finish_idx)),global_number_of_projected_points( &
                    & my_computational_node+1),mpi_integer,projected_face,global_number_of_projected_points, &
                    & global_mpi_displacements,mpi_integer,computational_environment%MPI_COMM,mpi_ierror) !PROJECTED_FACE
                  CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
                ENDIF
                DO ni=1,data_projection%NUMBER_OF_XI
                  CALL mpi_allgatherv(projected_xi(ni,sorting_ind_2(start_idx:finish_idx)),global_number_of_projected_points( &
                    & my_computational_node+1),mpi_double_precision,projected_xi(ni,:),global_number_of_projected_points, &
                    & global_mpi_displacements,mpi_double_precision,computational_environment%MPI_COMM,mpi_ierror) !PROJECTED_XI
                  CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
                ENDDO
                CALL mpi_allgatherv(projection_exit_tag(sorting_ind_2(start_idx:finish_idx)),global_number_of_projected_points( &
                  & my_computational_node+1),mpi_integer,projection_exit_tag,global_number_of_projected_points, &
                  & global_mpi_displacements,mpi_integer,computational_environment%MPI_COMM,mpi_ierror) !PROJECTION_EXIT_TAG
                CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
    DO ni=1,data_projection%COORDINATE_SYSTEM_DIMENSIONS
      CALL mpi_allgatherv(projection_vectors(ni, sorting_ind_2(start_idx:finish_idx)),global_number_of_projected_points( &
                    & my_computational_node+1),mpi_double_precision,projection_vectors(ni,:),global_number_of_projected_points, &
                    & global_mpi_displacements,mpi_double_precision,computational_environment%MPI_COMM,mpi_ierror)  !PROJECTION_VECTORS
                  CALL mpi_error_check("MPI_ALLGATHERV",mpi_ierror,err,error,*999)
    ENDDO
                !assign projection information to projected points
                DO data_point_idx=1,number_of_data_points
                  data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%EXIT_TAG=projection_exit_tag( &
                    & data_point_idx)
                  data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%ELEMENT_NUMBER=projected_element( &
                    & data_point_idx)
                  data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%DISTANCE=projected_distance(1,data_point_idx)
                  data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%XI(1:data_projection%NUMBER_OF_XI)= &
                    & projected_xi(1:data_projection%NUMBER_OF_XI,data_point_idx)
      data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%projectionVector( &
                    & 1:data_projection%COORDINATE_SYSTEM_DIMENSIONS)=projection_vectors( &
                    & 1:data_projection%COORDINATE_SYSTEM_DIMENSIONS,data_point_idx)
                ENDDO !data_point_idx
                projected_xi(:,sorting_ind_2)=projected_xi
    projection_vectors(:, sorting_ind_2)=projection_vectors
                projected_element(sorting_ind_2)=projected_element
                IF(data_projection%PROJECTION_TYPE==data_projection_boundary_lines_projection_type) THEN
                  DO data_point_idx=1,number_of_data_points
                    data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%ELEMENT_LINE_NUMBER=projected_face( &
                      & data_point_idx)
                  ENDDO !data_point_idx
                ELSEIF(data_projection%PROJECTION_TYPE==data_projection_boundary_faces_projection_type) THEN
                  DO data_point_idx=1,number_of_data_points
                    data_projection%DATA_PROJECTION_RESULTS(sorting_ind_2(data_point_idx))%ELEMENT_FACE_NUMBER=projected_face( &
                      & data_point_idx)
                  ENDDO !data_point_idx
                ENDIF
              ELSE !no need to use mpi
                SELECT CASE(data_projection%PROJECTION_TYPE)
                  CASE (data_projection_boundary_lines_projection_type) !Newton project to closest lines, and find miminum projection
                    DO data_point_idx=1,number_of_data_points
                      CALL data_projection_newton_lines_evaluate(data_projection,interpolated_point,data_points%DATA_POINTS( &
                        & data_point_idx)%position,closest_elements(data_point_idx,:),closest_faces(data_point_idx,:), &
                        & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%EXIT_TAG,data_projection% &
                        & data_projection_results(data_point_idx)%ELEMENT_NUMBER,data_projection% &
                        & data_projection_results(data_point_idx)%ELEMENT_LINE_NUMBER,data_projection%DATA_PROJECTION_RESULTS( &
                        & data_point_idx)%DISTANCE,data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI, &
                        & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector,err,error,*999)
                    ENDDO
                  CASE (data_projection_boundary_faces_projection_type) !find closest candidate faces
                    DO data_point_idx=1,number_of_data_points
                      CALL data_projection_newton_faces_evaluate(data_projection,interpolated_point,data_points%DATA_POINTS( &
                        & data_point_idx)%position,closest_elements(data_point_idx,:),closest_faces(data_point_idx,:), &
                        & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%EXIT_TAG,data_projection% &
                        & data_projection_results(data_point_idx)%ELEMENT_NUMBER,data_projection% &
                        & data_projection_results(data_point_idx)%ELEMENT_FACE_NUMBER,data_projection%DATA_PROJECTION_RESULTS( &
                        & data_point_idx)%DISTANCE,data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI, &
                        & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector,err,error,*999)
                    ENDDO
                  CASE (data_projection_all_elements_projection_type) !find closest candidate elements
                    SELECT CASE(data_projection%NUMBER_OF_XI)
                      CASE (1) !1D mesh
                        DO data_point_idx=1,number_of_data_points
                          CALL data_projection_newton_elements_evaluate_1(data_projection,interpolated_point,data_points% &
                            & data_points(data_point_idx)%position,closest_elements(data_point_idx,:),data_projection% &
                            & data_projection_results(data_point_idx)%EXIT_TAG,data_projection%DATA_PROJECTION_RESULTS( &
                            & data_point_idx)%ELEMENT_NUMBER, data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%DISTANCE, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector,&
                            & err,error,*999)
                        ENDDO
                      CASE (2) !2D mesh
                        DO data_point_idx=1,number_of_data_points
                          CALL data_projection_newton_elements_evaluate_2(data_projection,interpolated_point,data_points% &
                            & data_points(data_point_idx)%position,closest_elements(data_point_idx,:),data_projection% &
                            & data_projection_results(data_point_idx)%EXIT_TAG,data_projection%DATA_PROJECTION_RESULTS( &
                            & data_point_idx)%ELEMENT_NUMBER,data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%DISTANCE, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector, &
          & err,error,*999)
                        ENDDO
                      CASE (3) !3D mesh
                        DO data_point_idx=1,number_of_data_points
                          CALL data_projection_newton_elements_evaluate_3(data_projection,interpolated_point,data_points% &
                            & data_points(data_point_idx)%position,closest_elements(data_point_idx,:),data_projection% &
                            & data_projection_results(data_point_idx)%EXIT_TAG,data_projection%DATA_PROJECTION_RESULTS( &
                            & data_point_idx)%ELEMENT_NUMBER,data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%DISTANCE, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%XI, &
                            & data_projection%DATA_PROJECTION_RESULTS(data_point_idx)%projectionVector, &
          & err,error,*999)
                        ENDDO
                      CASE DEFAULT
                        CALL flag_error("Data projection number of xi is invalid",err,error,*999)
                    END SELECT !DATA_PROJECTION%NUMBER_OF_XI
                  CASE DEFAULT
                    CALL flag_error("No match for data projection type found",err,error,*999)
                END SELECT
              ENDIF !NUMBER_COMPUTATIONAL_NODES>1
              data_projection%DATA_PROJECTION_PROJECTED=.true.
            ELSE
              CALL flag_error("Data projection and projection field are not sharing the same mesh.",err,error,*999)
            ENDIF !ASSOCIATED(DATA_PROJECTION%MESH,PROJECTION_FIELD%DECOMPOSITION%MESH)
          ELSE
            CALL flag_error("Projection field have not been finished.",err,error,*999)
          ENDIF !PROJECTION_FIELD%FIELD_FINISHED
        ELSE
          CALL flag_error("Projection field is not associated.",err,error,*999)
        ENDIF !ASSOCIATED(PROJECTION_FIELD)
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF !DATA_PROJECTION%DATA_PROJECTION_FINISHED
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF !ASSOCIATED(DATA_PROJECTION)

    ! Deallocate arrays used within this routine
    IF(ALLOCATED(candidate_elements)) DEALLOCATE(candidate_elements)
    IF(ALLOCATED(candidate_faces)) DEALLOCATE(candidate_faces)
    IF(ALLOCATED(closest_elements)) DEALLOCATE(closest_elements)
    IF(ALLOCATED(closest_faces)) DEALLOCATE(closest_faces)
    IF(ALLOCATED(closest_distances)) DEALLOCATE(closest_distances)
    IF(ALLOCATED(global_to_local_number_of_closest_candidates)) DEALLOCATE(global_to_local_number_of_closest_candidates)
    IF(ALLOCATED(global_number_of_closest_candidates)) DEALLOCATE(global_number_of_closest_candidates)
    IF(ALLOCATED(global_mpi_displacements)) DEALLOCATE(global_mpi_displacements)
    IF(ALLOCATED(global_number_of_projected_points)) DEALLOCATE(global_number_of_projected_points)
    IF(ALLOCATED(projection_exit_tag)) DEALLOCATE(projection_exit_tag)
    IF(ALLOCATED(projected_element)) DEALLOCATE(projected_element)
    IF(ALLOCATED(projected_face)) DEALLOCATE(projected_face)
    IF(ALLOCATED(projected_distance)) DEALLOCATE(projected_distance)
    IF(ALLOCATED(projected_xi)) DEALLOCATE(projected_xi)
    IF(ALLOCATED(projection_vectors)) DEALLOCATE(projection_vectors)
    IF(ALLOCATED(sorting_ind_2)) DEALLOCATE(sorting_ind_2)
    IF(ALLOCATED(global_closest_distances)) DEALLOCATE(global_closest_distances)
    IF(ALLOCATED(sorting_ind_1)) DEALLOCATE(sorting_ind_1)

    exits("DataProjection_DataPointsProjectionEvaluate")
    RETURN
999 IF(ALLOCATED(candidate_elements)) DEALLOCATE(candidate_elements)
    IF(ALLOCATED(candidate_faces)) DEALLOCATE(candidate_faces)
    IF(ALLOCATED(closest_elements)) DEALLOCATE(closest_elements)
    IF(ALLOCATED(closest_faces)) DEALLOCATE(closest_faces)
    IF(ALLOCATED(closest_distances)) DEALLOCATE(closest_distances)
    IF(ALLOCATED(global_to_local_number_of_closest_candidates)) DEALLOCATE(global_to_local_number_of_closest_candidates)
    IF(ALLOCATED(global_number_of_closest_candidates)) DEALLOCATE(global_number_of_closest_candidates)
    IF(ALLOCATED(global_mpi_displacements)) DEALLOCATE(global_mpi_displacements)
    IF(ALLOCATED(global_number_of_projected_points)) DEALLOCATE(global_number_of_projected_points)
    IF(ALLOCATED(projection_exit_tag)) DEALLOCATE(projection_exit_tag)
    IF(ALLOCATED(projected_element)) DEALLOCATE(projected_element)
    IF(ALLOCATED(projected_face)) DEALLOCATE(projected_face)
    IF(ALLOCATED(projected_distance)) DEALLOCATE(projected_distance)
    IF(ALLOCATED(projected_xi)) DEALLOCATE(projected_xi)
    IF(ALLOCATED(projection_vectors)) DEALLOCATE(projection_vectors)
    IF(ALLOCATED(sorting_ind_2)) DEALLOCATE(sorting_ind_2)
    IF(ALLOCATED(global_closest_distances)) DEALLOCATE(global_closest_distances)
    IF(ALLOCATED(sorting_ind_1)) DEALLOCATE(sorting_ind_1)
998 errorsexits("DataProjection_DataPointsProjectionEvaluate",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_datapointsprojectionevaluate

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

  SUBROUTINE dataprojection_datapointspositionevaluate(dataProjection,field,fieldVariableType,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: dataProjection
    TYPE(field_type), POINTER :: field
    INTEGER(INTG), INTENT(IN) :: fieldVariableType
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: dataPointIdx,elementNumber,coordIdx
    TYPE(data_points_type), POINTER :: dataPoints
    TYPE(field_interpolated_point_ptr_type), POINTER :: interpolatedPoints(:)
    TYPE(field_interpolated_point_type), POINTER :: interpolatedPoint
    TYPE(field_interpolation_parameters_ptr_type), POINTER :: interpolationParameters(:)

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

    IF(ASSOCIATED(field)) THEN
      IF (field%TYPE==field_geometric_type.OR.field%TYPE==field_geometric_general_type) THEN
        IF(ASSOCIATED(dataprojection)) THEN
          datapoints=>dataprojection%DATA_POINTS
          IF(ASSOCIATED(datapoints)) THEN
            NULLIFY(interpolatedpoints)
            NULLIFY(interpolationparameters)
            CALL field_interpolation_parameters_initialise(field,interpolationparameters,err,error,*999, &
              & field_geometric_components_type)
            CALL field_interpolated_points_initialise(interpolationparameters,interpolatedpoints,err,error,*999, &
              & field_geometric_components_type)
            interpolatedpoint=>interpolatedpoints(fieldvariabletype)%PTR
            !Loop through data points
             DO datapointidx=1,datapoints%NUMBER_OF_DATA_POINTS
               elementnumber=dataprojection%DATA_PROJECTION_RESULTS(datapointidx)%ELEMENT_NUMBER
               CALL field_interpolation_parameters_element_get(field_values_set_type,elementnumber, &
                 & interpolationparameters(fieldvariabletype)%PTR,err,error,*999,field_geometric_components_type)
               CALL field_interpolate_xi(no_part_deriv,dataprojection%DATA_PROJECTION_RESULTS(datapointidx)%XI, &
                 & interpolatedpoint,err,error,*999,field_geometric_components_type)
               DO coordidx=1,SIZE(datapoints%DATA_POINTS(datapointidx)%position)
                 datapoints%DATA_POINTS(datapointidx)%position(coordidx)=interpolatedpoint%VALUES(coordidx,no_part_deriv)
               ENDDO !coordIdx
             ENDDO !dataPointIdx
           ELSE
             CALL flag_error("Data points is not associated.",err,error,*999)
           ENDIF
         ELSE
           CALL flag_error("Data projection is not associated.",err,error,*999)
         ENDIF
       ELSE
         CALL flag_error("Cannot evaluate data points position on field other than geometric or geometric general type.", &
           & err,error,*999)
       ENDIF
    ELSE
      CALL flag_error("Field is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_DataPointsPositionEvaluate")
    RETURN
999 errorsexits("DataProjection_DataPointsPositionEvaluate",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_datapointspositionevaluate

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

  SUBROUTINE dataprojection_maximuminterationupdateget(DATA_PROJECTION,MAXIMUM_ITERATION_UPDATE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(OUT) :: MAXIMUM_ITERATION_UPDATE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    enters("DataProjection_MaximumInterationUpdateGet",err,error,*999)

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        maximum_iteration_update=data_projection%MAXIMUM_ITERATION_UPDATE
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_MaximumInterationUpdateGet")
    RETURN
999 errorsexits("DataProjection_MaximumInterationUpdateGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_maximuminterationupdateget

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

  SUBROUTINE dataprojection_maximuminterationupdateset(DATA_PROJECTION,MAXIMUM_ITERATION_UPDATE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(IN) :: MAXIMUM_ITERATION_UPDATE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF((maximum_iteration_update>=0.1).AND.(maximum_iteration_update<=1)) THEN
          data_projection%MAXIMUM_ITERATION_UPDATE=maximum_iteration_update
        ELSE
          CALL flag_error("Data projection maximum iteration update must be between 0.1 and 1.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_MaximumInterationUpdateSet")
    RETURN
999 errorsexits("DataProjection_MaximumInterationUpdateSet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_maximuminterationupdateset


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

  SUBROUTINE dataprojection_maximumnumberofiterationsget(DATA_PROJECTION,MAXIMUM_NUMBER_OF_ITERATIONS,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: MAXIMUM_NUMBER_OF_ITERATIONS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        maximum_number_of_iterations=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_MaximumNumberOfIterationsGet")
    RETURN
999 errorsexits("DataProjection_MaximumNumberOfIterationsGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_maximumnumberofiterationsget

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

  SUBROUTINE dataprojection_maximumnumberofiterationsset(DATA_PROJECTION,MAXIMUM_NUMBER_OF_ITERATIONS,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: MAXIMUM_NUMBER_OF_ITERATIONS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF(maximum_number_of_iterations>=1) THEN
          data_projection%MAXIMUM_NUMBER_OF_ITERATIONS=maximum_number_of_iterations
        ELSE
          CALL flag_error("Data projection maximum number of iterations must be at least 1.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_MaximumNumberOfIterationsSet")
    RETURN
999 errorsexits("DataProjection_MaximumNumberOfIterationsSet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_maximumnumberofiterationsset

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

  SUBROUTINE data_projection_newton_elements_evaluate_1(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & projection_exit_tag,projection_element_number,projection_distance,projection_xi,projection_vector,err,error,*)
    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(1)
    REAL(DP), INTENT(OUT) :: PROJECTION_VECTOR(3)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    LOGICAL :: INSIDE_REGION,CONVERGED
    INTEGER(INTG) :: ELEMENT_NUMBER !local element number in current computational domain
    INTEGER(INTG) :: REGION_DIMENSIONS
    INTEGER(INTG) :: BOUND,EXIT_TAG
    REAL(DP) :: XI(1),XI_NEW(1),XI_UPDATE(1),XI_UPDATE_NORM
    REAL(DP) :: DISTANCE_VECTOR(3),RELATIVE_TOLERANCE,ABSOLUTE_TOLERANCE
    REAL(DP) :: FUNCTION_VALUE,FUNCTION_VALUE_NEW
    REAL(DP) :: FUNCTION_GRADIENT,FUNCTION_HESSIAN
    REAL(DP) :: MAXIMUM_DELTA,MINIMUM_DELTA,DELTA
    REAL(DP) :: PREDICTED_REDUCTION,PREDICTION_ACCURACY

    INTEGER(INTG) :: ne,itr1,itr2

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_exit_tag=data_projection_exit_tag_no_element
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
        maximum_delta=data_projection%MAXIMUM_ITERATION_UPDATE
        minimum_delta=0.025_dp*maximum_delta !need to set a minimum, in case if it gets too small
        DO ne=1,SIZE(candidate_elements,1) !project on each candidate elements
          element_number=candidate_elements(ne)
          exit_tag=data_projection_exit_tag_no_element
          converged=.false.
          delta=0.5_dp*maximum_delta !start at half the MAXIMUM_DELTA as we do not know if quadratic model is a good approximation yet
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number,interpolated_point% &
            & interpolation_parameters,err,error,*999,field_geometric_components_type)
          xi=data_projection%STARTING_XI
          CALL field_interpolate_xi(second_part_deriv,xi,interpolated_point,err,error,*999,field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values-interpolated_point%VALUES(:,no_part_deriv)
          function_value=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          main_loop: DO itr1=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(outer loop)
            !Check for bounds [0,1]
            IF(abs(xi(1))<zero_tolerance) THEN
              bound=-1 !bound at negative direction
            ELSEIF(abs(xi(1)-1.0_dp)<zero_tolerance) THEN
              bound=1 !bound at positive direction
            ELSE !inside the bounds
              bound=0
            ENDIF
            !FUNCTION_GRADIENT
            function_gradient=-2.0_dp* &
              & (dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,first_part_deriv)))
            !FUNCTION_HESSIAN
            function_hessian=-2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,second_part_deriv))- &
              & dot_product(interpolated_point%VALUES(:,first_part_deriv),interpolated_point%VALUES(:,first_part_deriv)))
            !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + EIGEN_SHIFT*I)g
            !The calculation of EIGEN_SHIFT are only approximated as oppose to the common trust region approach
            DO itr2=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations
              inside_region=.false.
              IF(function_hessian>absolute_tolerance) THEN !positive: minimum exists
                xi_update(1)=-function_gradient/function_hessian
                xi_update_norm=dabs(xi_update(1))
                inside_region=xi_update_norm<=delta
              ENDIF !positive
              IF(.NOT.inside_region) THEN !minimum not in the region
                xi_update(1)=-dsign(delta,function_gradient)
                xi_update_norm=delta
              ENDIF
              IF((bound/=0).AND.(bound>0.EQV.xi_update(1)>0.0_dp)) THEN !projection go out of element bound
                exit_tag=data_projection_exit_tag_bounds
                EXIT main_loop
              ENDIF
              converged=xi_update_norm<absolute_tolerance !first half of the convergence test (before collision detection)
              xi_new=xi+xi_update !update XI
              IF(xi_new(1)<0.0_dp) THEN !boundary collision check
                xi_new(1)=0.0_dp
              ELSEIF(xi_new(1)>1.0_dp) THEN
                xi_new(1)=1.0_dp
              ENDIF
              CALL field_interpolate_xi(second_part_deriv,xi_new,interpolated_point,err,error,*999,field_geometric_components_type)
              distance_vector=point_values-interpolated_point%VALUES(:,no_part_deriv)
              function_value_new=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
              converged=converged.AND.(dabs(function_value_new-function_value)/(1.0_dp+function_value)<relative_tolerance) !second half of the convergence test
              IF(converged) EXIT !converged: exit inner loop first
              IF((function_value_new-function_value)>absolute_tolerance) THEN !bad model: reduce step size
                IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                  exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                  EXIT main_loop
                ENDIF
                delta=dmax1(minimum_delta,0.25_dp*delta)
              ELSE
                predicted_reduction=xi_update(1)*(function_gradient+0.5_dp*function_hessian*xi_update(1))
                prediction_accuracy=(function_value_new-function_value)/predicted_reduction
                IF(prediction_accuracy<0.01_dp) THEN !bad model: reduce region size
                  IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                    exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                    EXIT main_loop
                  ENDIF
                  delta=dmax1(minimum_delta,0.5_dp*delta)
                ELSEIF(prediction_accuracy>0.9_dp.AND.prediction_accuracy<1.1_dp) THEN !good model: increase region size
                  delta=dmin1(maximum_delta,2.0_dp*delta)
                  EXIT
                ELSE !ok model: keep the current region size
                  EXIT
                ENDIF
              ENDIF
            ENDDO !itr2 (inner loop: adjust region size)
            function_value=function_value_new
            xi=xi_new
            IF(converged) THEN
              exit_tag=data_projection_exit_tag_converged
              EXIT
            ENDIF
          ENDDO main_loop !itr1 (outer loop)
          IF(exit_tag==data_projection_exit_tag_no_element.AND.itr1>=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS) &
            & exit_tag=data_projection_exit_tag_max_iteration
          IF((projection_exit_tag==data_projection_exit_tag_no_element).OR.(dsqrt(function_value)<projection_distance)) THEN
            projection_exit_tag=exit_tag
            projection_element_number=element_number
            projection_distance=dsqrt(function_value)
            projection_xi=xi
            projection_vector=distance_vector
          ENDIF
        ENDDO !ne
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_1")
    RETURN
999 errorsexits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_1",err,error)
    RETURN 1

  END SUBROUTINE data_projection_newton_elements_evaluate_1

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

  SUBROUTINE data_projection_newton_elements_evaluate_2(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & projection_exit_tag,projection_element_number,projection_distance,projection_xi,projection_vector,err,error,*)
    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(2)
    REAL(DP), INTENT(OUT) :: PROJECTION_VECTOR(3)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    TYPE(domain_mapping_type), POINTER :: DOMAIN_MAPPING
    LOGICAL :: FREE,CONVERGED,INSIDE_REGION
    INTEGER(INTG) :: ELEMENT_NUMBER
    INTEGER(INTG) :: MESH_COMPONENT_NUMBER,REGION_DIMENSIONS
    INTEGER(INTG) :: BOUND(2),EXIT_TAG
    REAL(DP) :: XI(2),XI_NEW(2),XI_UPDATE(2),XI_UPDATE_NORM
    REAL(DP) :: DISTANCE_VECTOR(3),RELATIVE_TOLERANCE,ABSOLUTE_TOLERANCE
    REAL(DP) :: FUNCTION_VALUE,FUNCTION_VALUE_NEW
    REAL(DP) :: FUNCTION_GRADIENT(2),FUNCTION_GRADIENT_NORM
    REAL(DP) :: FUNCTION_HESSIAN(2,2),HESSIAN_DIAGONAL(2)
    REAL(DP) :: TEMP1,TEMP2,DET,EIGEN_MIN,EIGEN_MAX,EIGEN_SHIFT
    REAL(DP) :: MAXIMUM_DELTA,MINIMUM_DELTA,DELTA
    REAL(DP) :: PREDICTED_REDUCTION,PREDICTION_ACCURACY


    INTEGER(INTG) :: ne,ni,nifix,itr1,itr2

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_exit_tag=data_projection_exit_tag_no_element
        mesh_component_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER
        domain_mapping=>interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%DOMAIN(mesh_component_number)% &
          & ptr%MAPPINGS%ELEMENTS
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
        maximum_delta=data_projection%MAXIMUM_ITERATION_UPDATE
        minimum_delta=0.025_dp*maximum_delta !need to set a minimum, in case if it gets too small
        DO ne=1,SIZE(candidate_elements,1) !project on each candidate elements
          element_number=candidate_elements(ne)
          exit_tag=data_projection_exit_tag_no_element
          converged=.false.
          delta=0.5_dp*maximum_delta !start at half the MAXIMUM_DELTA as we do not know if quadratic model is a good approximation yet
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          xi=data_projection%STARTING_XI
          CALL field_interpolate_xi(second_part_deriv,xi,interpolated_point,err,error,*999,field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values-interpolated_point%VALUES(:,no_part_deriv)
          function_value=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          main_loop: DO itr1=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(outer loop)
            !Check for bounds [0,1]
            DO ni=1,2
              IF(abs(xi(ni))<zero_tolerance) THEN
                bound(ni)=-1 !bound at negative direction
              ELSEIF(abs(xi(ni)-1.0_dp)<zero_tolerance) THEN
                bound(ni)=1 !bound at positive direction
              ELSE !inside the bounds
                bound(ni)=0
              ENDIF
            ENDDO !ni
            !FUNCTION_GRADIENT
            function_gradient(1)= &
              & -2.0_dp*(dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1)))
            function_gradient(2)= &
              & -2.0_dp*(dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s2)))
            !FUNCTION_HESSIAN
            function_hessian(1,1)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1_s1))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s1)))
            function_hessian(1,2)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s2)))
            function_hessian(2,2)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s2_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s2),interpolated_point%VALUES(:,part_deriv_s2)))
            !A model trust region approach, a Newton step is taken if the minimum lies inside the trust region (DELTA), if not, shift the step towards the steepest descent
            temp1=0.5_dp*(function_hessian(1,1)+function_hessian(2,2))
            temp2=dsqrt((0.5_dp*(function_hessian(1,1)-function_hessian(2,2)))**2+function_hessian(1,2)**2)
            eigen_min=temp1-temp2
            eigen_max=temp1+temp2
            function_gradient_norm=dsqrt(dot_product(function_gradient,function_gradient))
            DO itr2=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations
              temp1=function_gradient_norm/delta
              inside_region=(eigen_min>=temp1).AND.(eigen_min>absolute_tolerance) !estimate if the solution is inside the trust region without calculating a newton step, this also guarantees the hessian matrix is positive definite
              IF(inside_region) THEN
                det=eigen_min*eigen_max !det(H)
                hessian_diagonal(1)=function_hessian(1,1)
                hessian_diagonal(2)=function_hessian(2,2)
              ELSE
                eigen_shift=max(temp1-eigen_min,absolute_tolerance) !shift towards steepest decent
                det=temp1*(eigen_max+eigen_shift) !det(H)
                hessian_diagonal(1)=function_hessian(1,1)+eigen_shift
                hessian_diagonal(2)=function_hessian(2,2)+eigen_shift
              ENDIF
              xi_update(1)=-(hessian_diagonal(2)*function_gradient(1)-function_hessian(1,2)*function_gradient(2))/det
              xi_update(2)=(function_hessian(1,2)*function_gradient(1)-hessian_diagonal(1)*function_gradient(2))/det
              xi_update_norm=dsqrt(dot_product(xi_update,xi_update))
              free=.true.
              DO ni=1,2
                IF((bound(ni)/=0).AND.(bound(ni)>0.EQV.xi_update(ni)>0.0_dp)) THEN !projection go out of element bound
                  IF(.NOT.free) THEN !both xi are fixed
                    exit_tag=data_projection_exit_tag_bounds
                    EXIT main_loop
                  ENDIF
                  free=.false.
                  nifix=ni
                ENDIF
              ENDDO !ni
              IF(free) THEN !both xi are free
                IF(.NOT.inside_region) THEN
                  IF(xi_update_norm>0.0_dp) THEN
                    xi_update=delta/xi_update_norm*xi_update !readjust XI_UPDATE to lie on the region bound
                  ENDIF
                ENDIF
              ELSE !xi are not free
                xi_update(nifix)=0.0_dp
                ni=3-nifix
                inside_region=.false.
                IF(function_hessian(ni,ni)>0.0_dp) THEN !positive: minimum exists in the unbounded direction
                  xi_update(ni)=-function_gradient(ni)/function_hessian(ni,ni)
                  xi_update_norm=dabs(xi_update(ni))
                  inside_region=xi_update_norm<=delta
                ENDIF
                IF(.NOT.inside_region) THEN !minimum not in the region
                  xi_update(ni)=-dsign(delta,function_gradient(ni))
                  xi_update_norm=delta
                ENDIF
              ENDIF !if xi is free
              converged=xi_update_norm<absolute_tolerance !first half of the convergence test
              xi_new=xi+xi_update !update XI
              DO ni=1,2
                IF(xi_new(ni)<0.0_dp) THEN !boundary collision check
                  xi_new(ni)=0.0_dp
                  xi_new(3-ni)=xi(3-ni)-xi_update(3-ni)*xi(ni)/xi_update(ni)
                ELSEIF(xi_new(ni)>1.0_dp) THEN
                  xi_new(ni)=1.0_dp
                  xi_new(3-ni)=xi(3-ni)+xi_update(3-ni)*(1.0_dp-xi(ni))/xi_update(ni)
                ENDIF
              ENDDO
              CALL field_interpolate_xi(second_part_deriv,xi_new,interpolated_point,err,error,*999,field_geometric_components_type)
              distance_vector(1:region_dimensions)=point_values-interpolated_point%VALUES(:,no_part_deriv)
              function_value_new=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
              converged=converged.AND.(dabs(function_value_new-function_value)/(1.0_dp+function_value)<relative_tolerance) !second half of the convergence test (before collision detection)
              IF(converged) EXIT !converged: exit inner loop first
              IF((function_value_new-function_value)>absolute_tolerance) THEN !bad model: reduce step size
                IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                  exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                  EXIT main_loop
                ENDIF
                delta=dmax1(minimum_delta,0.25_dp*delta)
              ELSE
                predicted_reduction=dot_product(function_gradient,xi_update)+ &
                  & 0.5_dp*(xi_update(1)*(xi_update(1)*function_hessian(1,1)+2.0_dp*xi_update(2)*function_hessian(1,2))+ &
                  & xi_update(2)**2*function_hessian(2,2))
                prediction_accuracy=(function_value_new-function_value)/predicted_reduction
                IF(prediction_accuracy<0.01_dp) THEN !bad model: reduce region size
                  IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                    exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                    EXIT main_loop
                  ENDIF
                  delta=dmax1(minimum_delta,0.5_dp*delta)
                ELSEIF(prediction_accuracy>0.9_dp.AND.prediction_accuracy<1.1_dp) THEN !good model: increase region size
                  delta=dmin1(maximum_delta,2.0_dp*delta)
                  EXIT
                ELSE !ok model: keep the current region size
                  EXIT
                ENDIF
              ENDIF
            ENDDO !itr2 (inner loop: adjust region size)
            function_value=function_value_new
            xi=xi_new
            IF(converged) THEN
              exit_tag=data_projection_exit_tag_converged
              EXIT
            ENDIF
          ENDDO main_loop !itr1 (outer loop)
          IF(exit_tag==data_projection_exit_tag_no_element.AND.itr1>=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS) &
            & exit_tag=data_projection_exit_tag_max_iteration
          IF((projection_exit_tag==data_projection_exit_tag_no_element).OR.(dsqrt(function_value)<projection_distance)) THEN
            !IF(.NOT.ELEMENT_FOUND) ELEMENT_FOUND=.TRUE.
            projection_exit_tag=exit_tag
            projection_element_number=element_number
            projection_distance=dsqrt(function_value)
            projection_xi=xi
      projection_vector=distance_vector
          ENDIF
        ENDDO !ne
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_2")
    RETURN
999 errorsexits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_2",err,error)
    RETURN 1

  END SUBROUTINE data_projection_newton_elements_evaluate_2

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

  SUBROUTINE data_projection_newton_elements_evaluate_3(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & projection_exit_tag,projection_element_number,projection_distance,projection_xi,projection_vector,err,error,*)
    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(3)
    REAL(DP), INTENT(OUT) :: PROJECTION_VECTOR(3)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    TYPE(domain_mapping_type), POINTER :: DOMAIN_MAPPING
    LOGICAL :: FREE,CONVERGED,INSIDE_REGION
    INTEGER(INTG) :: ELEMENT_NUMBER
    INTEGER(INTG) :: MESH_COMPONENT_NUMBER
    INTEGER(INTG) :: NBOUND,BOUND(3),EXIT_TAG
    REAL(DP) :: XI(3),XI_NEW(3),XI_UPDATE(3),XI_UPDATE_NORM
    REAL(DP) :: DISTANCE_VECTOR(3),RELATIVE_TOLERANCE,ABSOLUTE_TOLERANCE
    REAL(DP) :: FUNCTION_VALUE,FUNCTION_VALUE_NEW
    REAL(DP) :: FUNCTION_GRADIENT(3),FUNCTION_GRADIENT_NORM,FUNCTION_GRADIENT2(2)
    REAL(DP) :: FUNCTION_HESSIAN(3,3),HESSIAN_DIAGONAL(3),FUNCTION_HESSIAN2(2,2),HESSIAN_DIAGONAL2(2)
    REAL(DP) :: TEMP1,TEMP2,TEMP3,TEMP4,DET,TRACE,TRACE2,EIGEN_MIN,EIGEN_MAX,EIGEN_SHIFT
    REAL(DP) :: MAXIMUM_DELTA,MINIMUM_DELTA,DELTA
    REAL(DP) :: PREDICTED_REDUCTION,PREDICTION_ACCURACY


    INTEGER(INTG) :: ne,ni,ni2(2),nb,nifix,nifix2(2),itr1,itr2,nbfix,mi

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_exit_tag=data_projection_exit_tag_no_element
        mesh_component_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER
        domain_mapping=>interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%DOMAIN(mesh_component_number)% &
          & ptr%MAPPINGS%ELEMENTS
        !REGION_DIMENSIONS=DATA_PROJECTION%COORDINATE_SYSTEM_DIMENSIONS
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
        maximum_delta=data_projection%MAXIMUM_ITERATION_UPDATE
        minimum_delta=0.025_dp*maximum_delta !need to set a minimum, in case if it gets too small
        DO ne=1,SIZE(candidate_elements,1) !project on each candidate elements
          element_number=candidate_elements(ne)
          exit_tag=data_projection_exit_tag_no_element
          converged=.false.
          delta=0.5_dp*maximum_delta !start at half the MAXIMUM_DELTA as we do not know if quadratic model is a good approximation yet
          CALL field_interpolation_parameters_element_get(field_values_set_type,element_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          xi=data_projection%STARTING_XI
          CALL field_interpolate_xi(second_part_deriv,xi,interpolated_point,err,error,*999,field_geometric_components_type)
          distance_vector=point_values-interpolated_point%VALUES(:,no_part_deriv)
          function_value=dot_product(distance_vector,distance_vector)
          main_loop: DO itr1=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(outer loop)
            !Check for bounds [0,1]
            DO ni=1,3
              IF(abs(xi(ni))<zero_tolerance) THEN
                bound(ni)=-1 !bound at negative direction
              ELSEIF(abs(xi(ni)-1.0_dp)<zero_tolerance) THEN
                bound(ni)=1 !bound at positive direction
              ELSE !inside the bounds
                bound(ni)=0
              ENDIF
            ENDDO !ni
            !FUNCTION_GRADIENT
            function_gradient(1)=-2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s1)))
            function_gradient(2)=-2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s2)))
            function_gradient(3)=-2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s3)))
            !FUNCTION_HESSIAN
            function_hessian(1,1)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s1_s1))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s1)))
            function_hessian(1,2)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s1_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s2)))
            function_hessian(1,3)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s1_s3))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s3)))
            function_hessian(2,2)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s2_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s2),interpolated_point%VALUES(:,part_deriv_s2)))
            function_hessian(2,3)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s2_s3))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s2),interpolated_point%VALUES(:,part_deriv_s3)))
            function_hessian(3,3)= -2.0_dp*(dot_product(distance_vector,interpolated_point%VALUES(:,part_deriv_s3_s3))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s3),interpolated_point%VALUES(:,part_deriv_s3)))
            !A model trust region approach, a Newton step is taken if the solution lies inside the trust region (DELTA), if not, shift the step towards the steepest descent
            trace=function_hessian(1,1)+function_hessian(2,2)+function_hessian(3,3) !tr(H)
            trace2=function_hessian(1,1)*function_hessian(2,2)+function_hessian(1,1)*function_hessian(3,3)+ &
              & function_hessian(2,2)*function_hessian(3,3)-function_hessian(1,2)**2-function_hessian(1,3)**2- &
              & function_hessian(2,3)**2 !tr(H**2)-(tr(H))**2
            det=function_hessian(1,1)*function_hessian(2,2)*function_hessian(3,3)- &
              & function_hessian(1,1)*function_hessian(2,3)**2-function_hessian(2,2)*function_hessian(1,3)**2- &
              & function_hessian(3,3)*function_hessian(1,2)**2+ &
              & 2.0_dp*function_hessian(1,2)*function_hessian(1,3)*function_hessian(2,3) !det(H)
            temp1=-trace/3.0_dp
            temp2=trace2/3.0_dp
            temp3=temp2-temp1**2 
            IF(temp3>-1.0e-5_dp) THEN !include some negatives for numerical errors
              eigen_min=-temp1 !all eigenvalues are the same
            ELSE
              temp3=dsqrt(-temp3)
              temp4=(det+3.0_dp*(temp1*temp2)-2.0_dp*temp1**3)/(2.0_dp*temp3**3)
              eigen_min=2.0_dp*temp3*dcos((dacos(temp4)+twopi)/3.0_dp)-temp1
            ENDIF
            function_gradient_norm=dsqrt(dot_product(function_gradient,function_gradient))
            DO itr2=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations
              temp1=function_gradient_norm/delta
              inside_region=(eigen_min>=temp1).AND.(eigen_min>absolute_tolerance) !estimate if the solution is inside the trust region without calculating a newton step, this also guarantees the hessian matrix is positive definite
              IF(inside_region) THEN
                hessian_diagonal(1)=function_hessian(1,1)
                hessian_diagonal(2)=function_hessian(2,2)
                hessian_diagonal(3)=function_hessian(3,3)
              ELSE
                eigen_shift=max(temp1-eigen_min,absolute_tolerance) !shift towards steepest decent
                det=det+eigen_shift*(trace2+eigen_shift*(trace+eigen_shift)) !shift the determinant
                hessian_diagonal(1)=function_hessian(1,1)+eigen_shift
                hessian_diagonal(2)=function_hessian(2,2)+eigen_shift
                hessian_diagonal(3)=function_hessian(3,3)+eigen_shift
              ENDIF
              temp2=function_hessian(1,3)*function_hessian(2,3)-function_hessian(1,2)*hessian_diagonal(3)
              temp3=function_hessian(1,2)*function_hessian(2,3)-function_hessian(1,3)*hessian_diagonal(2)
              temp4=function_hessian(1,2)*function_hessian(1,3)-function_hessian(2,3)*hessian_diagonal(1)
              xi_update(1)=((function_hessian(2,3)**2-hessian_diagonal(2)*hessian_diagonal(3))*function_gradient(1)- &
                & temp2*function_gradient(2)-temp3*function_gradient(3))/det
              xi_update(2)=((function_hessian(1,3)**2-hessian_diagonal(1)*hessian_diagonal(3))*function_gradient(2)- &
                & temp2*function_gradient(1)-temp4*function_gradient(3))/det
              xi_update(3)=((function_hessian(1,2)**2-hessian_diagonal(1)*hessian_diagonal(2))*function_gradient(3)- &
                & temp3*function_gradient(1)-temp4*function_gradient(2))/det
              xi_update_norm=dsqrt(dot_product(xi_update,xi_update))
              free=.true.
              nbound=0
              DO ni=1,3
                IF((bound(ni)/=0).AND.(bound(ni)>0.EQV.xi_update(ni)>0.0_dp)) THEN !projection go out of element bound
                  nbound=nbound+1
                  free=.false.
                  IF(nbound<=2) THEN
                    nifix2(nbound)=ni
                  ELSE !all xi are fixed
                    exit_tag=data_projection_exit_tag_bounds
                    EXIT main_loop
                  ENDIF
                ENDIF
              ENDDO !ni
              IF(free) THEN !all xi are free
                IF(.NOT.inside_region) THEN
                  IF(xi_update_norm>0.0_dp) THEN
                    xi_update=delta/xi_update_norm*xi_update !readjust XI_UPDATE to lie on the region bound
                  ENDIF
                ENDIF
              ELSE !at least one of the xi are not free
                !try 2D projection
                free=.true.
                nifix=nifix2(1)
                IF(nbound==2) THEN
                  IF(xi_update(nifix2(2))>xi_update(nifix2(1))) nifix=nifix2(2) !only fix the direction that is most strongly suggesting leaving the element
                ENDIF
                xi_update(nifix)=0.0_dp
                ni2(1)=1+mod(nifix,3)
                ni2(2)=1+mod(nifix+1,3)
                !FUNCTION_GRADIENT2
                function_gradient2(1)=function_gradient(ni2(1))
                function_gradient2(2)=function_gradient(ni2(2))
                !FUNCTION_HESSIAN2
                function_hessian2(1,1)=function_hessian(ni2(1),ni2(1))
                function_hessian2(1,2)=function_hessian(ni2(1),ni2(2))
                function_hessian2(2,2)=function_hessian(ni2(2),ni2(2))
                !re-estimate the trust solution in 2D
                temp1=0.5_dp*(function_hessian2(1,1)+function_hessian2(2,2))
                temp2=dsqrt((0.5_dp*(function_hessian2(1,1)-function_hessian2(2,2)))**2+function_hessian2(1,2)**2)
                eigen_min=temp1-temp2
                eigen_max=temp1+temp2
                temp3=dsqrt(dot_product(function_gradient2,function_gradient2))/delta
                inside_region=(eigen_min>=temp3).AND.(eigen_min>absolute_tolerance) !estimate if the solution is inside the trust region without calculating a newton step, this also guarantees the hessian matrix is positive definite
                IF(inside_region) THEN
                  det=eigen_min*eigen_max !determinant of FUNCTION_HESSIAN
                  hessian_diagonal2(1)=function_hessian2(1,1)
                  hessian_diagonal2(2)=function_hessian2(2,2)
                ELSE
                  eigen_shift=max(temp3-eigen_min,absolute_tolerance) !shift towards steepest decent
                  det=temp3*(eigen_max+eigen_shift) !determinant of shifted FUNCTION_HESSIAN
                  hessian_diagonal2(1)=function_hessian2(1,1)+eigen_shift
                  hessian_diagonal2(2)=function_hessian2(2,2)+eigen_shift
                ENDIF
                xi_update(ni2(1))=-(hessian_diagonal2(2)*function_gradient2(1)-function_hessian2(1,2)*function_gradient2(2))/det
                xi_update(ni2(2))=(function_hessian2(1,2)*function_gradient2(1)-hessian_diagonal2(1)*function_gradient2(2))/det
                xi_update_norm=dsqrt(dot_product(xi_update,xi_update))
                !check again for bounds
                DO nb=1,2
                  IF((bound(ni2(nb))/=0).AND.(bound(ni2(nb))>0.EQV.xi_update(ni2(nb))>0.0_dp)) THEN !projection go out of element bound
                    IF(.NOT.free) THEN !both xi are fixed
                      exit_tag=data_projection_exit_tag_bounds
                      EXIT main_loop
                    ENDIF
                    free=.false.
                    nbfix=nb
                  ENDIF
                ENDDO !ni
                IF(free) THEN !both xi are free
                  IF(.NOT.inside_region) THEN
                    IF(xi_update_norm>0.0_dp) THEN
                      xi_update=delta/xi_update_norm*xi_update !readjust XI_UPDATE to lie on the region bound
                    ENDIF
                  ENDIF
                ELSE !xi are not free
                  xi_update(ni2(nbfix))=0.0_dp
                  ni=ni2(3-nbfix)
                  inside_region=.false.
                  IF(function_hessian(ni,ni)>0.0_dp) THEN !positive: minimum exists in the unbounded direction
                    xi_update(ni)=-function_gradient(ni)/function_hessian(ni,ni)
                    xi_update_norm=dabs(xi_update(ni))
                    inside_region=xi_update_norm<=delta
                  ENDIF
                  IF(.NOT.inside_region) THEN !minimum not in the region
                    xi_update(ni)=-dsign(delta,function_gradient(ni))
                    xi_update_norm=delta
                  ENDIF
                ENDIF !if xi are free (2D)
              ENDIF !if xi are free (3D)
              converged=xi_update_norm<absolute_tolerance !first half of the convergence test
              xi_new=xi+xi_update !update XI
              DO ni=1,3
                IF(abs(xi_update(ni))<zero_tolerance) THEN !FPE Handling
                  IF(xi_new(ni)<0.0_dp) THEN
                    xi_new(ni)=0.0_dp
                    DO mi = 1,3
                      IF(mi /= ni) THEN
                        xi_new(mi)=xi(mi)-xi_update(mi)
                      ENDIF
                    ENDDO
                  ELSEIF(xi_new(ni)>1.0_dp) THEN
                    xi_new(ni)=1.0_dp
                    DO mi = 1,3
                      IF(mi /= ni) THEN
                        xi_new(mi)=xi(mi)+xi_update(mi)
                      ENDIF
                    ENDDO
                  ENDIF
                ELSEIF(xi_new(ni)<0.0_dp) THEN !boundary collision check
                  xi_new(ni)=0.0_dp
                  xi_new=xi-xi_update*xi(ni)/xi_update(ni)
                ELSEIF(xi_new(ni)>1.0_dp) THEN
                  xi_new(ni)=1.0_dp
                  xi_new=xi+xi_update*(1.0_dp-xi(ni))/xi_update(ni)
                ENDIF
              ENDDO !ni
              CALL field_interpolate_xi(second_part_deriv,xi_new,interpolated_point,err,error,*999,field_geometric_components_type)
              distance_vector=point_values-interpolated_point%VALUES(:,no_part_deriv)
              function_value_new=dot_product(distance_vector,distance_vector)
              converged=converged.AND.(dabs(function_value_new-function_value)/(1.0_dp+function_value)<relative_tolerance) !second half of the convergence test (before collision detection)
              IF(converged) EXIT !converged: exit inner loop first
              IF((function_value_new-function_value)>absolute_tolerance) THEN !bad model: reduce step size
                IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                  exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                  EXIT main_loop
                ENDIF
                delta=dmax1(minimum_delta,0.25_dp*delta)
              ELSE
                predicted_reduction=dot_product(function_gradient,xi_update)+ &
                  & 0.5_dp*(xi_update(1)*(xi_update(1)*function_hessian(1,1)+2.0_dp*xi_update(2)*function_hessian(1,2)+ &
                  & 2.0_dp*xi_update(3)*function_hessian(1,3))+xi_update(2)*(xi_update(2)*function_hessian(2,2)+ &
                  & 2.0_dp*xi_update(3)*function_hessian(2,3))+xi_update(2)**2*function_hessian(2,2))
                IF (abs(predicted_reduction) < zero_tolerance) THEN
                  converged = .true.  ! prediction reduction converged
                  EXIT
                ENDIF
                prediction_accuracy=(function_value_new-function_value)/predicted_reduction
                IF(prediction_accuracy<0.01_dp) THEN !bad model: reduce region size
                  IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                    exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                    EXIT main_loop
                  ENDIF
                  delta=dmax1(minimum_delta,0.5_dp*delta)
                ELSEIF(prediction_accuracy>0.9_dp.AND.prediction_accuracy<1.1_dp) THEN !good model: increase region size
                  delta=dmin1(maximum_delta,2.0_dp*delta)
                  EXIT
                ELSE !ok model: keep the current region size
                  EXIT
                ENDIF
              ENDIF
            ENDDO !itr2 (inner loop: adjust region size)
            function_value=function_value_new
            xi=xi_new
            IF(converged) THEN
              exit_tag=data_projection_exit_tag_converged
              EXIT
            ENDIF
          ENDDO main_loop !itr1 (outer loop)
          IF(exit_tag==data_projection_exit_tag_no_element.AND.itr1>=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS) &
            & exit_tag=data_projection_exit_tag_max_iteration
          IF((projection_exit_tag==data_projection_exit_tag_no_element).OR.(dsqrt(function_value)<projection_distance)) THEN
            !IF(.NOT.ELEMENT_FOUND) ELEMENT_FOUND=.TRUE.
            projection_exit_tag=exit_tag
            projection_element_number=element_number
            projection_distance=dsqrt(function_value)
            projection_xi=xi
      projection_vector=distance_vector
          ENDIF
        ENDDO !ne
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_3")
    RETURN
999 errorsexits("DATA_PROJECTION_NEWTON_ELEMENTS_EVALUATE_3",err,error)
    RETURN 1

  END SUBROUTINE data_projection_newton_elements_evaluate_3

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

  SUBROUTINE data_projection_newton_faces_evaluate(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & candidate_element_faces,projection_exit_tag,projection_element_number,projection_element_face_number,projection_distance, &
    & projection_xi,projection_vector,err,error,*)
    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENT_FACES(:)
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_FACE_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(2)
    REAL(DP), INTENT(OUT) :: PROJECTION_VECTOR(3)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    TYPE(domain_mapping_type), POINTER :: DOMAIN_MAPPING
    LOGICAL :: FREE,CONVERGED,INSIDE_REGION
    INTEGER(INTG) :: ELEMENT_NUMBER,ELEMENT_FACE_NUMBER,FACE_NUMBER
    INTEGER(INTG) :: MESH_COMPONENT_NUMBER,REGION_DIMENSIONS
    INTEGER(INTG) :: BOUND(2),EXIT_TAG
    REAL(DP) :: XI(2),XI_NEW(2),XI_UPDATE(2),XI_UPDATE_NORM
    REAL(DP) :: DISTANCE_VECTOR(3),RELATIVE_TOLERANCE,ABSOLUTE_TOLERANCE
    REAL(DP) :: FUNCTION_VALUE,FUNCTION_VALUE_NEW
    REAL(DP) :: FUNCTION_GRADIENT(2),FUNCTION_GRADIENT_NORM
    REAL(DP) :: FUNCTION_HESSIAN(2,2),HESSIAN_DIAGONAL(2)
    REAL(DP) :: TEMP1,TEMP2,DET,EIGEN_MIN,EIGEN_MAX,EIGEN_SHIFT
    REAL(DP) :: MAXIMUM_DELTA,MINIMUM_DELTA,DELTA
    REAL(DP) :: PREDICTED_REDUCTION,PREDICTION_ACCURACY


    INTEGER(INTG) :: ne,ni,nifix,itr1,itr2

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_exit_tag=data_projection_exit_tag_no_element
        mesh_component_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%MESH_COMPONENT_NUMBER
        domain_mapping=>interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%DOMAIN(mesh_component_number)% &
          & ptr%MAPPINGS%ELEMENTS
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
        maximum_delta=data_projection%MAXIMUM_ITERATION_UPDATE
        minimum_delta=0.025_dp*maximum_delta !need to set a minimum, in case if it gets too small
        DO ne=1,SIZE(candidate_elements,1) !project on each candidate elements
          element_number=candidate_elements(ne)
          element_face_number=candidate_element_faces(ne)
          face_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_FACES(element_face_number)
          exit_tag=data_projection_exit_tag_no_element
          converged=.false.
          delta=0.5_dp*maximum_delta !start at half the MAXIMUM_DELTA as we do not know if quadratic model is a good approximation yet
          CALL field_interpolation_parameters_face_get(field_values_set_type,face_number, &
            & interpolated_point%INTERPOLATION_PARAMETERS,err,error,*999,field_geometric_components_type)
          xi=data_projection%STARTING_XI
          CALL field_interpolate_xi(second_part_deriv,xi,interpolated_point,err,error,*999,field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values-interpolated_point%VALUES(:,no_part_deriv)
          function_value=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          main_loop: DO itr1=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(outer loop)
            !Check for bounds [0,1]
            DO ni=1,2
              IF(abs(xi(ni))<zero_tolerance) THEN
                bound(ni)=-1 !bound at negative direction
              ELSEIF(abs(xi(ni)-1.0_dp)<zero_tolerance) THEN
                bound(ni)=1 !bound at positive direction
              ELSE !inside the bounds
                bound(ni)=0
              ENDIF
            ENDDO !ni
            !FUNCTION_GRADIENT
            function_gradient(1)= &
              & -2.0_dp*(dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1)))
            function_gradient(2)= &
              & -2.0_dp*(dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s2)))
            !FUNCTION_HESSIAN
            function_hessian(1,1)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1_s1))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s1)))
            function_hessian(1,2)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s1_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s1),interpolated_point%VALUES(:,part_deriv_s2)))
            function_hessian(2,2)= -2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,part_deriv_s2_s2))- &
              & dot_product(interpolated_point%VALUES(:,part_deriv_s2),interpolated_point%VALUES(:,part_deriv_s2)))
            !A model trust region approach, a Newton step is taken if the minimum lies inside the trust region (DELTA), if not, shift the step towards the steepest descent
            temp1=0.5_dp*(function_hessian(1,1)+function_hessian(2,2))
            temp2=dsqrt((0.5_dp*(function_hessian(1,1)-function_hessian(2,2)))**2+function_hessian(1,2)**2)
            eigen_min=temp1-temp2
            eigen_max=temp1+temp2
            function_gradient_norm=dsqrt(dot_product(function_gradient,function_gradient))
            DO itr2=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations
              temp1=function_gradient_norm/delta
              inside_region=(eigen_min>=temp1).AND.(eigen_min>absolute_tolerance) !estimate if the solution is inside the trust region without calculating a newton step, this also guarantees the hessian matrix is positive definite
              IF(inside_region) THEN
                det=eigen_min*eigen_max !det(H)
                hessian_diagonal(1)=function_hessian(1,1)
                hessian_diagonal(2)=function_hessian(2,2)
              ELSE
                eigen_shift=max(temp1-eigen_min,absolute_tolerance) !shift towards steepest decent
                det=temp1*(eigen_max+eigen_shift) !det(H)
                hessian_diagonal(1)=function_hessian(1,1)+eigen_shift
                hessian_diagonal(2)=function_hessian(2,2)+eigen_shift
              ENDIF
              xi_update(1)=-(hessian_diagonal(2)*function_gradient(1)-function_hessian(1,2)*function_gradient(2))/det
              xi_update(2)=(function_hessian(1,2)*function_gradient(1)-hessian_diagonal(1)*function_gradient(2))/det
              xi_update_norm=dsqrt(dot_product(xi_update,xi_update))
              free=.true.
              DO ni=1,2
                IF((bound(ni)/=0).AND.(bound(ni)>0.EQV.xi_update(ni)>0.0_dp)) THEN !projection go out of element bound
                  IF(.NOT.free) THEN !both xi are fixed
                    exit_tag=data_projection_exit_tag_bounds
                    EXIT main_loop
                  ENDIF
                  free=.false.
                  nifix=ni
                ENDIF
              ENDDO !ni
              IF(free) THEN !both xi are free
                IF(.NOT.inside_region) THEN
                  IF(xi_update_norm>0.0_dp) THEN
                    xi_update=delta/xi_update_norm*xi_update !readjust XI_UPDATE to lie on the region bound
                  ENDIF
                ENDIF
              ELSE !xi are not free
                xi_update(nifix)=0.0_dp
                ni=3-nifix
                inside_region=.false.
                IF(function_hessian(ni,ni)>0.0_dp) THEN !positive: minimum exists in the unbounded direction
                  xi_update(ni)=-function_gradient(ni)/function_hessian(ni,ni)
                  xi_update_norm=dabs(xi_update(ni))
                  inside_region=xi_update_norm<=delta
                ENDIF
                IF(.NOT.inside_region) THEN !minimum not in the region
                  xi_update(ni)=-dsign(delta,function_gradient(ni))
                  xi_update_norm=delta
                ENDIF
              ENDIF !if xi is free
              converged=xi_update_norm<absolute_tolerance !first half of the convergence test
              xi_new=xi+xi_update !update XI
              DO ni=1,2
                IF(xi_new(ni)<0.0_dp) THEN !boundary collision check
                  xi_new(ni)=0.0_dp
                  xi_new(3-ni)=xi(3-ni)-xi_update(3-ni)*xi(ni)/xi_update(ni)
                ELSEIF(xi_new(ni)>1.0_dp) THEN
                  xi_new(ni)=1.0_dp
                  xi_new(3-ni)=xi(3-ni)+xi_update(3-ni)*(1.0_dp-xi(ni))/xi_update(ni)
                ENDIF
              ENDDO
              CALL field_interpolate_xi(second_part_deriv,xi_new,interpolated_point,err,error,*999,field_geometric_components_type)
              distance_vector=point_values-interpolated_point%VALUES(:,no_part_deriv)
              function_value_new=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
              converged=converged.AND.(dabs(function_value_new-function_value)/(1.0_dp+function_value)<relative_tolerance) !second half of the convergence test (before collision detection)
              IF(converged) THEN
                function_value=function_value_new
                EXIT !converged: exit inner loop first
              ENDIF
              IF((function_value_new-function_value)>absolute_tolerance) THEN !bad model: reduce step size
                IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                  exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                  function_value=function_value_new
                  EXIT main_loop
                ENDIF
                delta=dmax1(minimum_delta,0.25_dp*delta)
              ELSE
                predicted_reduction=dot_product(function_gradient,xi_update)+ &
                  & 0.5_dp*(xi_update(1)*(xi_update(1)*function_hessian(1,1)+2.0_dp*xi_update(2)*function_hessian(1,2))+ &
                  & xi_update(2)**2*function_hessian(2,2))
                prediction_accuracy=(function_value_new-function_value)/predicted_reduction
                IF(prediction_accuracy<0.01_dp) THEN !bad model: reduce region size
                  IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                    exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                    function_value=function_value_new
                    EXIT main_loop
                  ENDIF
                  delta=dmax1(minimum_delta,0.5_dp*delta)
                ELSEIF(prediction_accuracy>0.9_dp.AND.prediction_accuracy<1.1_dp) THEN !good model: increase region size
                  delta=dmin1(maximum_delta,2.0_dp*delta)
                  function_value=function_value_new
                  EXIT
                ELSE !ok model: keep the current region size
                  function_value=function_value_new
                  EXIT
                ENDIF
              ENDIF
            ENDDO !itr2 (inner loop: adjust region size)
            function_value=function_value_new
            xi=xi_new
            IF(converged) THEN
              exit_tag=data_projection_exit_tag_converged
              EXIT
            ENDIF
          ENDDO main_loop !itr1 (outer loop)
          IF(exit_tag==data_projection_exit_tag_no_element.AND.itr1>=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS) &
            & exit_tag=data_projection_exit_tag_max_iteration
          IF((projection_exit_tag==data_projection_exit_tag_no_element).OR.(dsqrt(function_value)<projection_distance)) THEN
            !IF(.NOT.ELEMENT_FOUND) ELEMENT_FOUND=.TRUE.
            projection_exit_tag=exit_tag
            projection_element_number=element_number
            projection_element_face_number=element_face_number
            projection_distance=dsqrt(function_value)
            projection_xi=xi
      projection_vector=distance_vector
          ENDIF
        ENDDO !ne
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_NEWTON_FACES_EVALUATE")
    RETURN
999 errorsexits("DATA_PROJECTION_NEWTON_FACES_EVALUATE",err,error)
    RETURN 1

  END SUBROUTINE data_projection_newton_faces_evaluate

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

  SUBROUTINE data_projection_newton_lines_evaluate(DATA_PROJECTION,INTERPOLATED_POINT,POINT_VALUES,CANDIDATE_ELEMENTS, &
    & candidate_element_lines,projection_exit_tag,projection_element_number,projection_element_line_number,projection_distance, &
    & projection_xi,projection_vector,err,error,*)
    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
    REAL(DP), INTENT(IN) :: POINT_VALUES(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENTS(:)
    INTEGER(INTG), INTENT(IN) :: CANDIDATE_ELEMENT_LINES(:)
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_LINE_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(1)
    REAL(DP), INTENT(OUT) :: PROJECTION_VECTOR(3)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR

    !Local Variables
    LOGICAL :: INSIDE_REGION,CONVERGED
    INTEGER(INTG) :: ELEMENT_NUMBER,ELEMENT_LINE_NUMBER,LINE_NUMBER !local element number in current computational domain
    INTEGER(INTG) :: REGION_DIMENSIONS
    INTEGER(INTG) :: BOUND,EXIT_TAG
    REAL(DP) :: XI(1),XI_NEW(1),XI_UPDATE(1),XI_UPDATE_NORM
    REAL(DP) :: DISTANCE_VECTOR(3),RELATIVE_TOLERANCE,ABSOLUTE_TOLERANCE
    REAL(DP) :: FUNCTION_VALUE,FUNCTION_VALUE_NEW
    REAL(DP) :: FUNCTION_GRADIENT,FUNCTION_HESSIAN
    REAL(DP) :: MAXIMUM_DELTA,MINIMUM_DELTA,DELTA
    REAL(DP) :: PREDICTED_REDUCTION,PREDICTION_ACCURACY

    INTEGER(INTG) :: ne,itr1,itr2

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_exit_tag=data_projection_exit_tag_no_element
        region_dimensions=data_projection%COORDINATE_SYSTEM_DIMENSIONS
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
        absolute_tolerance=data_projection%ABSOLUTE_TOLERANCE
        maximum_delta=data_projection%MAXIMUM_ITERATION_UPDATE
        minimum_delta=0.025_dp*maximum_delta !need to set a minimum, in case if it gets too small
        DO ne=1,SIZE(candidate_elements,1) !project on each candidate elements
          element_number=candidate_elements(ne)
          element_line_number=candidate_element_lines(ne)
          line_number=interpolated_point%INTERPOLATION_PARAMETERS%FIELD%DECOMPOSITION%TOPOLOGY%ELEMENTS%ELEMENTS( &
            & element_number)%ELEMENT_LINES(element_line_number)
          exit_tag=data_projection_exit_tag_no_element
          converged=.false.
          delta=0.5_dp*maximum_delta !start at half the MAXIMUM_DELTA as we do not know if quadratic model is a good approximation yet
          CALL field_interpolation_parameters_line_get(field_values_set_type,line_number,interpolated_point% &
            & interpolation_parameters,err,error,*999,field_geometric_components_type)
          xi=data_projection%STARTING_XI
          CALL field_interpolate_xi(second_part_deriv,xi,interpolated_point,err,error,*999,field_geometric_components_type)
          distance_vector(1:region_dimensions)=point_values-interpolated_point%VALUES(:,no_part_deriv)
          function_value=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
          main_loop: DO itr1=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(outer loop)
            !Check for bounds [0,1]
            IF(abs(xi(1))<zero_tolerance) THEN
              bound=-1 !bound at negative direction
            ELSEIF(abs(xi(1)-1.0_dp)<zero_tolerance) THEN
              bound=1 !bound at positive direction
            ELSE !inside the bounds
              bound=0
            ENDIF
            !FUNCTION_GRADIENT
            function_gradient=-2.0_dp* &
              & (dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,first_part_deriv)))
            !FUNCTION_HESSIAN
            function_hessian=-2.0_dp*(&
              & dot_product(distance_vector(1:region_dimensions),interpolated_point%VALUES(:,second_part_deriv))- &
              & dot_product(interpolated_point%VALUES(:,first_part_deriv),interpolated_point%VALUES(:,first_part_deriv)))
            !A model trust region approach, directly taken from CMISS CLOS22: V = -(H + EIGEN_SHIFT*I)g
            !The calculation of EIGEN_SHIFT are only approximated as oppose to the common trust region approach
            DO itr2=1,data_projection%MAXIMUM_NUMBER_OF_ITERATIONS !(inner loop: adjust region size) usually EXIT at 1 or 2 iterations
              inside_region=.false.
              IF(function_hessian>absolute_tolerance) THEN !positive: minimum exists
                xi_update(1)=-function_gradient/function_hessian
                xi_update_norm=dabs(xi_update(1))
                inside_region=xi_update_norm<=delta
              ENDIF !positive
              IF(.NOT.inside_region) THEN !minimum not in the region
                xi_update(1)=-dsign(delta,function_gradient)
                xi_update_norm=delta
              ENDIF
              IF((bound/=0).AND.(bound>0.EQV.xi_update(1)>0.0_dp)) THEN !projection go out of element bound
                exit_tag=data_projection_exit_tag_bounds
                EXIT main_loop
              ENDIF
              converged=xi_update_norm<absolute_tolerance !first half of the convergence test (before collision detection)
              xi_new=xi+xi_update !update XI
              IF(xi_new(1)<0.0_dp) THEN !boundary collision check
                xi_new(1)=0.0_dp
              ELSEIF(xi_new(1)>1.0_dp) THEN
                xi_new(1)=1.0_dp
              ENDIF
              CALL field_interpolate_xi(second_part_deriv,xi_new,interpolated_point,err,error,*999,field_geometric_components_type)
              distance_vector=point_values-interpolated_point%VALUES(:,no_part_deriv)
              function_value_new=dot_product(distance_vector(1:region_dimensions),distance_vector(1:region_dimensions))
              converged=converged.AND.(dabs(function_value_new-function_value)/(1.0_dp+function_value)<relative_tolerance) !second half of the convergence test
              IF(converged) EXIT !converged: exit inner loop first
              IF((function_value_new-function_value)>absolute_tolerance) THEN !bad model: reduce step size
                IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                  exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                  EXIT main_loop
                ENDIF
                delta=dmax1(minimum_delta,0.25_dp*delta)
              ELSE
                predicted_reduction=xi_update(1)*(function_gradient+0.5_dp*function_hessian*xi_update(1))
                prediction_accuracy=(function_value_new-function_value)/predicted_reduction
                IF(prediction_accuracy<0.01_dp) THEN !bad model: reduce region size
                  IF(delta<=minimum_delta) THEN !something went wrong, MINIMUM_DELTA too large? not likely to happen if MINIMUM_DELTA is small
                    exit_tag=data_projection_exit_tag_max_iteration ! it will get stucked!!
                    EXIT main_loop
                  ENDIF
                  delta=dmax1(minimum_delta,0.5_dp*delta)
                ELSEIF(prediction_accuracy>0.9_dp.AND.prediction_accuracy<1.1_dp) THEN !good model: increase region size
                  delta=dmin1(maximum_delta,2.0_dp*delta)
                  EXIT
                ELSE !ok model: keep the current region size
                  EXIT
                ENDIF
              ENDIF
            ENDDO !itr2 (inner loop: adjust region size)
            function_value=function_value_new
            xi=xi_new
            IF(converged) THEN
              exit_tag=data_projection_exit_tag_converged
              EXIT
            ENDIF
          ENDDO main_loop !itr1 (outer loop)
          IF(exit_tag==data_projection_exit_tag_no_element.AND.itr1>=data_projection%MAXIMUM_NUMBER_OF_ITERATIONS) &
            & exit_tag=data_projection_exit_tag_max_iteration
          IF((projection_exit_tag==data_projection_exit_tag_no_element).OR.(dsqrt(function_value)<projection_distance)) THEN
            projection_exit_tag=exit_tag
            projection_element_number=element_number
            projection_element_line_number=element_line_number
            projection_distance=dsqrt(function_value)
            projection_xi=xi
      projection_vector=distance_vector
          ENDIF
        ENDDO !ne
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_NEWTON_LINES_EVALUATE")
    RETURN
999 errorsexits("DATA_PROJECTION_NEWTON_LINES_EVALUATE",err,error)
    RETURN 1

  END SUBROUTINE data_projection_newton_lines_evaluate

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

  SUBROUTINE dataprojection_numberofclosestelementsget(DATA_PROJECTION,NUMBER_OF_CLOSEST_ELEMENTS,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: NUMBER_OF_CLOSEST_ELEMENTS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        number_of_closest_elements=data_projection%NUMBER_OF_CLOSEST_ELEMENTS
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_NumberOfClosestElementsGet")
    RETURN
999 errorsexits("DataProjection_NumberOfClosestElementsGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_numberofclosestelementsget

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

  SUBROUTINE dataprojection_numberofclosestelementsset(DATA_PROJECTION,NUMBER_OF_CLOSEST_ELEMENTS,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: NUMBER_OF_CLOSEST_ELEMENTS
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF(number_of_closest_elements>=1) THEN
          data_projection%NUMBER_OF_CLOSEST_ELEMENTS=number_of_closest_elements
        ELSE
          CALL flag_error("Data projection number of closest elements must be at least 1.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_NumberOfClosestElementsSet")
    RETURN
999 errorsexits("DataProjection_NumberOfClosestElementsSet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_numberofclosestelementsset

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

  SUBROUTINE dataprojection_projectioncandidatesset(dataProjection,elementUserNumber,localFaceLineNumbers,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: dataProjection
    INTEGER(INTG), INTENT(IN) :: elementUserNumber(:)
    INTEGER(INTG), INTENT(IN) :: localFaceLineNumbers(:)
    INTEGER(INTG), INTENT(OUT) :: err
    TYPE(varying_string), INTENT(OUT) :: error
    !Local Variables
    INTEGER(INTG) :: elementIdx,elementMeshNumber
    INTEGER(INTG) :: meshComponentNumber=1
    LOGICAL :: elementExists

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

    IF(ASSOCIATED(dataprojection)) THEN
      IF(SIZE(elementusernumber,1)==SIZE(localfacelinenumbers,1)) THEN
        ALLOCATE(dataprojection%candidateElementNumbers(SIZE(elementusernumber,1)),stat=err)
        IF(err/=0) CALL flag_error("Could not allocate candidiate element numbers.",err,error,*998)
        ALLOCATE(dataprojection%localFaceLineNumbers(SIZE(localfacelinenumbers,1)),stat=err)
        IF(err/=0) CALL flag_error("Could not allocate candidiate local face/line numbers.",err,error,*999)
        DO elementidx=1,SIZE(elementusernumber,1)
          CALL meshtopologyelementcheckexists(dataprojection%MESH,meshcomponentnumber,elementusernumber(elementidx), &
            & elementexists,elementmeshnumber,err,error,*999)
          IF(elementexists) THEN
            dataprojection%candidateElementNumbers(elementidx)=elementusernumber(elementidx)
            dataprojection%localFaceLineNumbers(elementidx)=localfacelinenumbers(elementidx)
          ELSE
            CALL flag_error("Element with user number ("//trim(number_to_vstring &
              & (elementusernumber(elementidx),"*",err,error))//") does not exist.",err,error,*999)
          ENDIF
        ENDDO !elementIdx
      ELSE
        CALL flag_error("Input user element numbers and face numbers sizes do not match.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_ProjectionCandidatesSet")
    RETURN
999 IF(ALLOCATED(dataprojection%candidateElementNumbers)) THEN
      DEALLOCATE(dataprojection%candidateElementNumbers)
    END IF
    IF(ALLOCATED(dataprojection%localFaceLineNumbers)) THEN
      DEALLOCATE(dataprojection%localFaceLineNumbers)
    END IF
998 errorsexits("DataProjection_ProjectionCandidatesSet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_projectioncandidatesset

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

  SUBROUTINE data_projection_projection_type_get(DATA_PROJECTION,PROJECTION_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        projection_type=data_projection%PROJECTION_TYPE
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_PROJECTION_TYPE_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_PROJECTION_TYPE_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_projection_type_get

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

  SUBROUTINE data_projection_projection_type_set(DATA_PROJECTION,PROJECTION_TYPE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: PROJECTION_TYPE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    REAL(DP), ALLOCATABLE :: STARTING_XI(:)

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        data_projection%PROJECTION_TYPE=projection_type
        SELECT CASE(projection_type)
          CASE (data_projection_boundary_lines_projection_type)
            data_projection%NUMBER_OF_XI=1
          CASE (data_projection_boundary_faces_projection_type)
            data_projection%NUMBER_OF_XI=2
          CASE (data_projection_all_elements_projection_type)
            data_projection%NUMBER_OF_XI=data_projection%MESH%NUMBER_OF_DIMENSIONS
          CASE DEFAULT
            CALL flag_error("Input projection type is undefined.",err,error,*999)
        END SELECT
        IF(data_projection%NUMBER_OF_XI/=SIZE(data_projection%STARTING_XI,1)) THEN
          ALLOCATE(starting_xi(data_projection%NUMBER_OF_XI),stat=err)
          IF(err/=0) CALL flag_error("Could not allocate starting xi.",err,error,*999)
          IF(data_projection%NUMBER_OF_XI>SIZE(data_projection%STARTING_XI,1)) THEN
            starting_xi(1:SIZE(data_projection%STARTING_XI,1))=data_projection%STARTING_XI(1:SIZE(data_projection%STARTING_XI,1))
            starting_xi(SIZE(data_projection%STARTING_XI,1):data_projection%NUMBER_OF_XI)=0.5_dp
          ELSE
            starting_xi(1:SIZE(data_projection%STARTING_XI,1))=data_projection%STARTING_XI(1:data_projection%NUMBER_OF_XI)
          ENDIF
          DEALLOCATE(data_projection%STARTING_XI)
          ALLOCATE(data_projection%STARTING_XI(data_projection%NUMBER_OF_XI),stat=err)
          IF(err/=0) CALL flag_error("Could not allocate data projection starting xi.",err,error,*999)
          data_projection%STARTING_XI(1:data_projection%NUMBER_OF_XI)=starting_xi(1:data_projection%NUMBER_OF_XI)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_PROJECTION_TYPE_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_PROJECTION_TYPE_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_projection_type_set

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

  SUBROUTINE data_projection_relative_tolerance_get(DATA_PROJECTION,RELATIVE_TOLERANCE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(OUT) :: RELATIVE_TOLERANCE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        relative_tolerance=data_projection%RELATIVE_TOLERANCE
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RELATIVE_TOLERANCE_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_RELATIVE_TOLERANCE_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_relative_tolerance_get

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

  SUBROUTINE data_projection_relative_tolerance_set(DATA_PROJECTION,RELATIVE_TOLERANCE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(IN) :: RELATIVE_TOLERANCE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    enters("DATA_PROJECTION_RELATIVE_TOLERANCE_SET",err,error,*999)

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF(relative_tolerance>=0) THEN
          data_projection%RELATIVE_TOLERANCE=relative_tolerance
        ELSE
          CALL flag_error("Data projection relative tolerance must be a positive real number.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RELATIVE_TOLERANCE_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_RELATIVE_TOLERANCE_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_relative_tolerance_set


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

  SUBROUTINE data_projection_starting_xi_get(DATA_PROJECTION,STARTING_XI,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(OUT) :: STARTING_XI(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    CHARACTER(LEN=MAXSTRLEN) :: LOCAL_ERROR

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(SIZE(starting_xi,1)>=SIZE(data_projection%STARTING_XI,1)) THEN
          starting_xi(1:SIZE(data_projection%STARTING_XI,1))=data_projection%STARTING_XI(1:SIZE(data_projection%STARTING_XI,1))
        ELSE
          WRITE(local_error,'("The size of the supplied starting xi  array of ",I2," is too small. The size must be >= ",I2,".")' &
            & )SIZE(starting_xi,1),SIZE(data_projection%STARTING_XI,1)
          CALL flag_error(local_error,err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_STARTING_XI_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_STARTING_XI_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_starting_xi_get

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

  SUBROUTINE data_projection_starting_xi_set(DATA_PROJECTION,STARTING_XI,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    REAL(DP), INTENT(IN) :: STARTING_XI(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: ni
    LOGICAL :: VALID_INPUT

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        CALL flag_error("Data projection have been finished.",err,error,*999)
      ELSE
        IF(SIZE(starting_xi,1)==SIZE(data_projection%STARTING_XI,1)) THEN
          valid_input=.true.
          DO ni=1,SIZE(starting_xi,1)
            IF((starting_xi(ni)<0).OR.(starting_xi(ni)>1)) THEN
              valid_input=.false.
              EXIT !this do
            ENDIF
          ENDDO
          IF(valid_input) THEN
            data_projection%STARTING_XI(1:SIZE(starting_xi))=starting_xi(1:SIZE(starting_xi))
          ELSE
            CALL flag_error("Data projection starting xi must be between 0 and 1.",err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection starting xi dimension mismatch.",err,error,*999)
        ENDIF
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_STARTING_XI_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_STARTING_XI_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_starting_xi_set

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

  SUBROUTINE data_projection_element_set(DATA_PROJECTION,DATA_POINT_USER_NUMBER,ELEMENT_USER_NUMBER,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    INTEGER(INTG), INTENT(IN) :: ELEMENT_USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER
    LOGICAL :: DATA_POINT_EXISTS

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

    IF(ASSOCIATED(data_projection)) THEN
      CALL dataprojection_datapointcheckexist(data_projection,data_point_user_number,data_point_exists, &
        & data_point_global_number,err,error,*999)
      IF(data_point_exists) THEN
        data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%ELEMENT_NUMBER=element_user_number
      ELSE
        CALL flag_error("Data point with user number ("//trim(number_to_vstring &
            & (data_point_user_number,"*",err,error))//") does not exist.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_ELEMENT_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_ELEMENT_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_element_set

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

  SUBROUTINE data_projection_result_distance_get(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_DISTANCE,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_DISTANCE
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          projection_distance=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%DISTANCE
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RESULT_DISTANCE_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_RESULT_DISTANCE_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_result_distance_get

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

  SUBROUTINE data_projection_label_get_vs(DATA_PROJECTION,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(varying_string), INTENT(OUT) :: LABEL
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      label=data_projection%LABEL
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_LABEL_GET_VS")
    RETURN
999 errorsexits("DATA_PROJECTION_LABEL_GET_VS",err,error)
    RETURN 1

  END SUBROUTINE data_projection_label_get_vs

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

  SUBROUTINE data_projection_label_get_c(DATA_PROJECTION,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    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("DATA_PROJECTION_LABEL_GET_C",err,error,*999)

    IF(ASSOCIATED(data_projection)) THEN
      c_length=len(label)
      vs_length=len_trim(data_projection%LABEL)
      IF(c_length>vs_length) THEN
        label=char(len_trim(data_projection%LABEL))
      ELSE
        label=char(data_projection%LABEL,c_length)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_LABEL_GET_C")
    RETURN
999 errorsexits("DATA_PROJECTION_LABEL_GET_C",err,error)
    RETURN 1

  END SUBROUTINE data_projection_label_get_c

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

  SUBROUTINE data_projection_label_set_c(DATA_PROJECTION,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    CHARACTER(LEN=*), INTENT(IN) :: LABEL
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      data_projection%LABEL=label
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_LABEL_SET_C")
    RETURN
999 errorsexits("DATA_PROJECTION_LABEL_SET_C",err,error)
    RETURN 1

  END SUBROUTINE data_projection_label_set_c

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

  SUBROUTINE data_projection_label_set_vs(DATA_PROJECTION,LABEL,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    TYPE(varying_string), INTENT(IN) :: LABEL
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables

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

    IF(ASSOCIATED(data_projection)) THEN
      data_projection%LABEL=label
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_LABEL_SET_VS")
    RETURN
999 errorsexits("DATA_PROJECTION_LABEL_SET_VS",err,error)
    RETURN 1

  END SUBROUTINE data_projection_label_set_vs

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

  SUBROUTINE data_projection_result_element_number_get(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_ELEMENT_NUMBER, &
    & err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          projection_element_number=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%ELEMENT_NUMBER
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RESULT_ELEMENT_NUMBER_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_RESULT_ELEMENT_NUMBER_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_result_element_number_get

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

  SUBROUTINE dataprojection_resultelementfacenumberget(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_ELEMENT_FACE_NUMBER &
    & ,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_FACE_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          ! Check if boundary faces projection type was set
          IF(data_projection%PROJECTION_TYPE==data_projection_boundary_faces_projection_type) THEN
            CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
              & data_point_global_number,err,error,*999)
            projection_element_face_number=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%ELEMENT_FACE_NUMBER
          ELSE
            CALL flag_error("Data projection data projection projection type is not set to boundary faces projection type.", &
              & err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_ResultElementFaceNumberGet")
    RETURN
999 errorsexits("DataProjection_ResultElementFaceNumberGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_resultelementfacenumberget

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

  SUBROUTINE dataprojection_resultelementlinenumberget(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_ELEMENT_LINE_NUMBER &
    & ,err,error,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_ELEMENT_LINE_NUMBER
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          ! Check if boundary lines projection type was set
          IF(data_projection%PROJECTION_TYPE==data_projection_boundary_lines_projection_type) THEN
            CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
              & data_point_global_number,err,error,*999)
            projection_element_line_number=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%ELEMENT_LINE_NUMBER
          ELSE
            CALL flag_error("Data projection data projection projection type is not set to boundary lines projection type.", &
              & err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_ResultElementLineNumberGet")
    RETURN
999 errorsexits("DataProjection_ResultElementLineNumberGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_resultelementlinenumberget

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

  SUBROUTINE data_projection_result_exit_tag_get(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_EXIT_TAG,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    INTEGER(INTG), INTENT(OUT) :: PROJECTION_EXIT_TAG
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          projection_exit_tag=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%EXIT_TAG
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RESULT_EXIT_TAG_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_RESULT_EXIT_TAG_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_result_exit_tag_get

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

  SUBROUTINE data_projection_result_xi_get(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_XI,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    REAL(DP), INTENT(OUT) :: PROJECTION_XI(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          IF(SIZE(projection_xi,1)==SIZE(data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%XI,1)) THEN
            projection_xi=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%XI
          ELSE
            CALL flag_error("projection xi has size of "//trim(number_to_vstring(SIZE(projection_xi,1),"*",err,error))// &
              & "but it needs to have size of "// &
              & trim(number_to_vstring(SIZE(data_projection%DATA_PROJECTION_RESULTS &
              & (data_point_global_number)%XI,1),"*",err,error))// &
              & "." ,err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RESULT_XI_GET")
    RETURN
999 errorsexits("DATA_PROJECTION_RESULT_XI_GET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_result_xi_get

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

  SUBROUTINE data_projection_result_xi_set(DATA_PROJECTION,DATA_POINT_USER_NUMBER,PROJECTION_XI,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    REAL(DP), INTENT(IN) :: PROJECTION_XI(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          IF(SIZE(projection_xi,1)==SIZE(data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%XI,1)) THEN
            data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%XI(1:SIZE(projection_xi,1))= &
              & projection_xi(1:SIZE(projection_xi,1))
          ELSE
            CALL flag_error("projection xi has size of "//trim(number_to_vstring(SIZE(projection_xi,1),"*",err,error))// &
              & "but it needs to have size of "// &
              & trim(number_to_vstring(SIZE(data_projection%DATA_PROJECTION_RESULTS &
              & (data_point_global_number)%XI,1),"*",err,error))// &
              & "." ,err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DATA_PROJECTION_RESULT_XI_SET")
    RETURN
999 errorsexits("DATA_PROJECTION_RESULT_XI_SET",err,error)
    RETURN 1

  END SUBROUTINE data_projection_result_xi_set

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

  SUBROUTINE dataprojection_resultprojectionvectorget(DATA_PROJECTION,DATA_POINT_USER_NUMBER,projectionVector,ERR,ERROR,*)

    !Argument variables
    TYPE(data_projection_type), POINTER :: DATA_PROJECTION
    INTEGER(INTG), INTENT(IN) :: DATA_POINT_USER_NUMBER
    REAL(DP), INTENT(OUT) :: projectionVector(:)
    INTEGER(INTG), INTENT(OUT) :: ERR
    TYPE(varying_string), INTENT(OUT) :: ERROR
    !Local Variables
    INTEGER(INTG) :: DATA_POINT_GLOBAL_NUMBER

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

    IF(ASSOCIATED(data_projection)) THEN
      IF(data_projection%DATA_PROJECTION_FINISHED) THEN
        IF(data_projection%DATA_PROJECTION_PROJECTED) THEN
          CALL dataprojection_datapointglobalnumberget(data_projection,data_point_user_number, &
            & data_point_global_number,err,error,*999)
          IF(SIZE(projectionvector,1)>=SIZE(data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)%projectionVector, &
            & 1)) THEN
            projectionvector=data_projection%DATA_PROJECTION_RESULTS(data_point_global_number)% &
              & projectionvector(1:data_projection%COORDINATE_SYSTEM_DIMENSIONS)
          ELSE
            CALL flag_error("projection vector has size of "//trim(number_to_vstring(SIZE(projectionvector,1),"*",err,error))// &
              & "but it needs to have size of "// &
              & trim(number_to_vstring(SIZE(data_projection%DATA_PROJECTION_RESULTS &
              & (data_point_global_number)%projectionVector,1),"*",err,error))// &
              & "." ,err,error,*999)
          ENDIF
        ELSE
          CALL flag_error("Data projection have not been projected.",err,error,*999)
        ENDIF
      ELSE
        CALL flag_error("Data projection have not been finished.",err,error,*999)
      ENDIF
    ELSE
      CALL flag_error("Data projection is not associated.",err,error,*999)
    ENDIF

    exits("DataProjection_ResultProjectionVectorGet")
    RETURN
999 errorsexits("DataProjection_ResultProjectionVectorGet",err,error)
    RETURN 1

  END SUBROUTINE dataprojection_resultprojectionvectorget

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

END MODULE data_projection_routines