59 type(
field_t),
intent(inout) :: field_data
60 class(*),
intent(in) :: object
61 real(kind=
dp),
intent(in),
optional :: max_distance
63 real(kind=
dp) :: max_dist
65 if (
present(max_distance))
then
66 max_dist = max_distance
68 max_dist = huge(0.0_dp)
76 call neko_error(
"signed_distance_field: Object type not supported.")
95 type(
field_t),
intent(inout) :: field_data
97 real(kind=
dp),
intent(in) :: max_distance
102 real(kind=
dp),
dimension(3) :: p
103 real(kind=
dp) :: distance
106 field_data%x = 0.0_dp
107 total_size = field_data%dof%size()
109 call search_tree%init(
mesh%nelv)
110 call search_tree%build(
mesh%el)
112 if (search_tree%get_size() .ne.
mesh%nelv)
then
113 call neko_error(
"signed_distance_field_tri_mesh: &
114 & Error building the search tree.")
117 do id = 1, total_size
118 p(1) = field_data%dof%x(id, 1, 1, 1)
119 p(2) = field_data%dof%y(id, 1, 1, 1)
120 p(3) = field_data%dof%z(id, 1, 1, 1)
124 field_data%x(id, 1, 1, 1) =
real(distance, kind=
rp)
147 real(kind=
dp),
intent(in) :: p(3)
148 real(kind=
dp),
intent(in) :: max_distance
151 real(kind=
dp) :: distance, weighted_sign
152 real(kind=
dp) :: cd, cs
153 real(kind=
dp) :: tol = 1e-6_dp
156 weighted_sign = 0.0_dp
162 if (abs(cd - distance) / distance .lt. tol)
then
163 weighted_sign = weighted_sign + cs
164 else if (cd .lt. distance)
then
168 distance = min(cd, distance)
171 distance = sign(min(distance, max_distance), weighted_sign)
195 class(
tri_t),
dimension(:),
intent(in) :: object_list
196 real(kind=
dp),
dimension(3),
intent(in) :: p
197 real(kind=
dp),
intent(in) :: max_distance
199 real(kind=
dp) :: distance
200 real(kind=
dp) :: weighted_sign
202 real(kind=
dp),
parameter :: tol = 1.0e-6_dp
205 integer :: current_index
208 type(
aabb_t) :: current_aabb
209 integer :: current_object_index
210 real(kind=
dp) :: current_distance
211 real(kind=
dp) :: current_sign
217 type(
aabb_t) :: search_box
219 integer :: root_index, left_index, right_index
220 real(kind=
dp) :: random_value
223 call simple_stack%init(
size(object_list) * 2)
224 call search_box%init(p - max_distance, p + max_distance)
228 root_index = tree%get_root_index()
229 root_box = tree%get_aabb(root_index)
231 if (.not. root_box%overlaps(search_box))
then
232 distance = max_distance
233 weighted_sign = 1.0_dp
238 call random_number(random_value)
239 current_object_index = floor(random_value *
size(object_list) + 1)
240 call element_distance(object_list(current_object_index), p, distance, weighted_sign)
241 distance = distance + object_list(current_object_index)%diameter()
244 call search_box%init(p - distance, p + distance)
245 call simple_stack%push(root_index)
248 do while (.not. simple_stack%is_empty())
249 current_index = simple_stack%pop()
252 current_node = tree%get_node(current_index)
253 current_aabb = current_node%get_aabb()
255 if (current_node%is_leaf())
then
256 if (distance .lt. current_node%min_distance(p))
then
260 current_object_index = current_node%get_object_index()
261 call element_distance(object_list(current_object_index), p, current_distance, current_sign)
264 if (abs(current_distance - distance) / distance .lt. tol)
then
265 weighted_sign = weighted_sign + current_sign
266 else if (current_distance .lt. distance)
then
267 weighted_sign = current_sign
270 distance = min(distance, current_distance)
273 if (distance .gt. current_aabb%get_diameter())
then
274 call search_box%init(p - distance, p + distance)
278 left_index = tree%get_left_index(current_index)
280 left_node = tree%get_left_node(current_index)
281 if (left_node%aabb%overlaps(search_box))
then
282 call simple_stack%push(left_index)
286 right_index = tree%get_right_index(current_index)
288 right_node = tree%get_right_node(current_index)
289 if (right_node%aabb%overlaps(search_box))
then
290 call simple_stack%push(right_index)
296 if (distance .gt. max_distance)
then
297 distance = max_distance
299 distance = sign(distance, weighted_sign)
309 class(*),
intent(in) ::
element
310 real(kind=
dp),
dimension(3),
intent(in) :: p
311 real(kind=
dp),
intent(out) :: distance
312 real(kind=
dp),
intent(out),
optional :: weighted_sign
319 print *,
"Error: Element type not supported."
345 type(
tri_t),
intent(in) :: triangle
346 real(kind=
dp),
dimension(3),
intent(in) :: p
348 real(kind=
dp),
intent(out) :: distance
349 real(kind=
dp),
intent(out),
optional :: weighted_sign
351 real(kind=
dp),
dimension(3) :: v1, v2, v3
352 real(kind=
dp),
dimension(3) :: normal
353 real(kind=
dp) :: normal_length
354 real(kind=
dp) :: b1, b2, b3
357 real(kind=
dp),
dimension(3) :: edge
358 real(kind=
dp) :: tol = 1e-10_dp
360 real(kind=
dp) :: face_distance
364 v1 = triangle%pts(1)%p%x
365 v2 = triangle%pts(2)%p%x
366 v3 = triangle%pts(3)%p%x
368 normal =
cross(v2 - v1, v3 - v1)
369 normal_length = norm2(normal)
371 if (normal_length .lt. tol)
then
372 distance = huge(1.0_dp)
373 weighted_sign = 0.0_dp
376 normal = normal / normal_length
380 face_distance = dot_product(p - v1, normal)
383 b1 = dot_product(normal,
cross(v2 - v1,
projection - v1)) / normal_length
384 b2 = dot_product(normal,
cross(v3 - v2,
projection - v2)) / normal_length
385 b3 = dot_product(normal,
cross(v1 - v3,
projection - v3)) / normal_length
387 if (b1 .le. tol)
then
389 t = dot_product(p - v1, edge) / norm2(edge)
390 t =
max(0.0_dp, min(1.0_dp, t))
393 else if (b2 .le. tol)
then
395 t = dot_product(p - v2, edge) / norm2(edge)
396 t =
max(0.0_dp, min(1.0_dp, t))
399 else if (b3 .le. tol)
then
401 t = dot_product(p - v3, edge) / norm2(edge)
402 t =
max(0.0_dp, min(1.0_dp, t))
408 if (
present(weighted_sign))
then
409 weighted_sign = face_distance / distance
411 distance = sign(distance, face_distance)
421 real(kind=
dp),
dimension(3),
intent(in) :: a
422 real(kind=
dp),
dimension(3),
intent(in) :: b
423 real(kind=
dp),
dimension(3) :: c
425 c(1) = a(2) * b(3) - a(3) * b(2)
426 c(2) = a(3) * b(1) - a(1) * b(3)
427 c(3) = a(1) * b(2) - a(2) * b(1)
Axis Aligned Bounding Box (aabb) Tree data structure.
integer, parameter, public aabb_null_node
Axis Aligned Bounding Box (aabb) implementation in Fortran.
integer, parameter, public dp
integer, parameter, public rp
Global precision used in computations.
Project x onto X, the space of old solutions and back again.
Module containing Signed Distance Functions.
subroutine signed_distance_field_tri_mesh(field_data, mesh, max_distance)
Signed distance field for a triangular mesh.
pure real(kind=dp) function, dimension(3) cross(a, b)
Compute cross product of two vectors.
subroutine element_distance_triangle(triangle, p, distance, weighted_sign)
Signed distance function for a triangle.
subroutine signed_distance_field(field_data, object, max_distance)
Signed distance field.
real(kind=dp) function tri_mesh_aabb_tree(tree, object_list, p, max_distance)
Signed distance function using an AABB tree.
subroutine element_distance(element, p, distance, weighted_sign)
Main interface for the signed distance function for an element.
real(kind=dp) function tri_mesh_brute_force(mesh, p, max_distance)
Signed distance function.
Implements a dynamic stack ADT.
Defines a triangular surface mesh.
Defines a triangular element.
Axis Aligned Bounding Box (aabb) data structure.
Node type for the Axis Aligned Bounding Box (aabb) Tree.
Axis Aligned Bounding Box (aabb) Tree.
A point in with coordinates .