@@ -967,7 +967,7 @@ contains
967967 real (wp), intent (out ) :: distance
968968
969969 integer :: i
970- real (wp) :: dist_min, dist_proj, dist_v1, dist_v2 , t
970+ real (wp) :: dist_min, dist , t
971971 real (wp) :: v1(1 :2 ), v2(1 :2 ), edge(1 :2 ), pv(1 :2 )
972972 real (wp) :: edge_len_sq, proj(1 :2 )
973973
@@ -982,11 +982,9 @@ contains
982982 v2(2 ) = boundary_v(i, 2 , 2 , pid)
983983
984984 ! Edge vector and point- to - v1 vector
985- edge(1 ) = v2(1 ) - v1(1 )
986- edge(2 ) = v2(2 ) - v1(2 )
985+ edge = v2 - v1
987986 pv(1 ) = point(1 ) - v1(1 )
988987 pv(2 ) = point(2 ) - v1(2 )
989-
990988 edge_len_sq = dot_product (edge, edge)
991989
992990 ! Parameter of projection onto the edge line
@@ -998,31 +996,18 @@ contains
998996
999997 ! Check if projection falls within the segment
1000998 if (t >= 0._wp .and. t <= 1._wp ) then
1001- proj(1 ) = v1(1 ) + t* edge(1 )
1002- proj(2 ) = v1(2 ) + t* edge(2 )
1003- dist_proj = sqrt ((point(1 ) - proj(1 ))** 2 + (point(2 ) - proj(2 ))** 2 )
1004-
1005- if (dist_proj < dist_min) then
1006- dist_min = dist_proj
1007- normals(1 ) = boundary_v(i, 3 , 1 , pid)
1008- normals(2 ) = boundary_v(i, 3 , 2 , pid)
1009- end if
1010- else
1011- ! Check distance to first vertex
1012- dist_v1 = sqrt (pv(1 )* pv(1 ) + pv(2 )* pv(2 ))
1013- if (dist_v1 < dist_min) then
1014- dist_min = dist_v1
1015- normals(1 ) = boundary_v(i, 3 , 1 , pid)
1016- normals(2 ) = boundary_v(i, 3 , 2 , pid)
1017- end if
999+ proj = v1 + t* edge
1000+ dist = sqrt ((point(1 ) - proj(1 ))** 2 + (point(2 ) - proj(2 ))** 2 )
1001+ else if (t < 0._wp ) then ! negative t means that v1 is the closest point on the edge
1002+ dist = sqrt ((point(1 ) - v1(1 ))** 2 + (point(2 ) - v1(2 ))** 2 )
1003+ else ! t > 1 means that v2 is the closest point on the line edge
1004+ dist = sqrt ((point(1 ) - v2(1 ))** 2 + (point(2 ) - v2(2 ))** 2 )
1005+ end if
10181006
1019- ! Check distance to second vertex
1020- dist_v2 = sqrt ((point(1 ) - v2(1 ))** 2 + (point(2 ) - v2(2 ))** 2 )
1021- if (dist_v2 < dist_min) then
1022- dist_min = dist_v2
1023- normals(1 ) = boundary_v(i, 3 , 1 , pid)
1024- normals(2 ) = boundary_v(i, 3 , 2 , pid)
1025- end if
1007+ if (dist < dist_min) then
1008+ dist_min = dist
1009+ normals(1 ) = boundary_v(i, 3 , 1 , pid)
1010+ normals(2 ) = boundary_v(i, 3 , 2 , pid)
10261011 end if
10271012 end do
10281013
0 commit comments