From 7899e5a04628c012fd7f15ada522f22641e27695 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 15 May 2026 09:15:31 +0100 Subject: [PATCH 1/4] Use a bisection algorithm to try to speed up `get_ielement()` for larger values of nelement. Leave original implementation in comments. --- src/Interpolation.jl | 79 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 64 insertions(+), 15 deletions(-) diff --git a/src/Interpolation.jl b/src/Interpolation.jl index 2437055..eba1e7e 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -1,4 +1,4 @@ -export value_in_coordinate_domain, +export value_in_coordinate_domain, get_ielement, interpolate_1D @@ -15,7 +15,7 @@ function value_in_coordinate_domain(xval::Float64, xcoord::FiniteElementCoordina # boundary points elseif (abs(xval-xebs[1]) < 100*tolerance) || (abs(xval-xebs[end]) < 100*tolerance) in_domain = true - else + else in_domain = false end return in_domain @@ -25,27 +25,76 @@ end Function to find the element in which the value xval sits in the assembled grid of the coordinate xcoord. """ +# linear search algorithm +# function get_ielement(xval::Float64, +# xcoord::FiniteElementCoordinate) +# @inbounds begin +# xebs = xcoord.element_boundaries +# nelement = xcoord.nelement +# tolerance = 1.0e-14 +# ielement = -1 +# # find the element +# for j in 1:nelement +# # check for internal points +# if (xval - xebs[j])*(xebs[j+1] - xval) > tolerance +# ielement = j +# break +# # check for boundary points +# elseif (abs(xval-xebs[j]) < 100*tolerance) || (abs(xval-xebs[j+1]) < 100*tolerance && j == nelement) +# ielement = j +# break +# end +# end +# if ielement < 1 +# error("xval=$xval is not within the coordinate $(xcoord.name)") +# end +# return ielement +# end +# end function get_ielement(xval::Float64, xcoord::FiniteElementCoordinate) + # search by bisection @inbounds begin xebs = xcoord.element_boundaries nelement = xcoord.nelement - tolerance = 1.0e-14 - ielement = -1 - # find the element - for j in 1:nelement - # check for internal points - if (xval - xebs[j])*(xebs[j+1] - xval) > tolerance - ielement = j - break - # check for boundary points - elseif (abs(xval-xebs[j]) < 100*tolerance) || (abs(xval-xebs[j+1]) < 100*tolerance && j == nelement) + tolerance = eps(Float64) + i = 1 # lower limit + k = nelement+1 # upper limit + # check extreme limits before main bisection loop + xi = xebs[i] - xval + xk = xebs[k] - xval + if abs(xi) < tolerance + # xval is xebs[i], so in element i + return i + elseif abs(xk) < tolerance + # xval is xebs[k], so in element k - 1 + return k - 1 # ielement range is 1:nelement, but we have nelement+1 boundaries + elseif xk*xi > tolerance # positive if root not in [xebs[i],xebs[k]] + error("xval=$xval is not within the coordinate $(xcoord.name)") + end + j = i + div(k - i, 2) # midpoint + ielement = i + # find the element by bisection + # note that because the element range i 1:nelement, but we have + # nelement+1 element boundaries, we default to returning the + # lower limit index in the loops below + while k - i > 1 + xj = xebs[j] - xval + if abs(xj) < tolerance + # ielement = j - div(j , nelement + 1) + # j is always < nelement + 1, so no need for additional div(j, nelement + 1) term ielement = j break + elseif xj*(xval - xebs[k]) > tolerance + # root between j and k + i = j # set lower limit index to j + j = i + div(k - i, 2) # set new midpoint < k + else + # root between i and j + k = j # set upper limit to j + j = i + div(k - i, 2) # set new midpoint =< i end - end - if ielement < 1 - error("xval=$xval is not within the coordinate $(xcoord.name)") + ielement = i end return ielement end From e203ad780a3f38eb3d1953e69f510422791e4759 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 12 Jun 2026 14:15:49 +0100 Subject: [PATCH 2/4] Tighten tolerances in `value_in_coordinate_domain()`. --- src/Interpolation.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Interpolation.jl b/src/Interpolation.jl index eba1e7e..923ecbc 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -8,12 +8,12 @@ the domain covered by the coordinate grid `xcoord`. """ function value_in_coordinate_domain(xval::Float64, xcoord::FiniteElementCoordinate) xebs = xcoord.element_boundaries - tolerance = 1.0e-14 + tolerance = eps(Float64) # internal point if (xval - xebs[1])*(xebs[end] - xval) > tolerance in_domain = true # boundary points - elseif (abs(xval-xebs[1]) < 100*tolerance) || (abs(xval-xebs[end]) < 100*tolerance) + elseif (abs(xval-xebs[1]) < tolerance) || (abs(xval-xebs[end]) < tolerance) in_domain = true else in_domain = false From 5532579755c21f73e008c5707ab24d3b0c2fb755 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 15 Jun 2026 09:01:36 +0100 Subject: [PATCH 3/4] Check lower and upper internal endpoints of the reference grids are -1 and 1, respectively, to make sure that the continuous Galerkin method can use an assembly for the global grid and operators of the type implemented. --- src/FiniteElementAssembly.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/FiniteElementAssembly.jl b/src/FiniteElementAssembly.jl index ebabd4d..f35e088 100644 --- a/src/FiniteElementAssembly.jl +++ b/src/FiniteElementAssembly.jl @@ -364,11 +364,20 @@ struct FiniteElementCoordinate{Tbc <: AbstractBoundaryCondition} nelement = length(element_data) # number of grid points per element ngrid = length(element_data[1].lpoly_data.x_nodes) - # check ngrid the same for each element + # checks on inputs for j in 2:nelement + # check ngrid the same for each element if !(length(element_data[j].lpoly_data.x_nodes) == ngrid) error("length(element_data[j].lpoly_data.x_nodes) /= ngrid \n Number of nodes in reference grid must be the same for each element") end + # check lower internal endpoint = -1 + if !(abs(element_data[j].lpoly_data.x_nodes[1] + 1.0) < eps(Float64)) + error("Internal element boundaries must be included in the domain with element_data[$j].lpoly_data.x_nodes[1] = -1.0") + end + # check upper internal endpoint = 1 + if !(abs(element_data[j-1].lpoly_data.x_nodes[end] - 1.0) < eps(Float64)) + error("Internal element boundaries must be included in the domain with element_data[$(j-1)].lpoly_data.x_nodes[end] = 1.0") + end end # extract shift, scale, and boundary values element_scale = Array{Float64,1}(undef,nelement) From f8f2a2f0d75ca5ee089555667052b04bbe51380f Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 17 Jun 2026 08:23:14 +0100 Subject: [PATCH 4/4] Removed commented `get_ielement()` implementation. --- src/Interpolation.jl | 26 -------------------------- 1 file changed, 26 deletions(-) diff --git a/src/Interpolation.jl b/src/Interpolation.jl index 923ecbc..e10bb1c 100644 --- a/src/Interpolation.jl +++ b/src/Interpolation.jl @@ -25,32 +25,6 @@ end Function to find the element in which the value xval sits in the assembled grid of the coordinate xcoord. """ -# linear search algorithm -# function get_ielement(xval::Float64, -# xcoord::FiniteElementCoordinate) -# @inbounds begin -# xebs = xcoord.element_boundaries -# nelement = xcoord.nelement -# tolerance = 1.0e-14 -# ielement = -1 -# # find the element -# for j in 1:nelement -# # check for internal points -# if (xval - xebs[j])*(xebs[j+1] - xval) > tolerance -# ielement = j -# break -# # check for boundary points -# elseif (abs(xval-xebs[j]) < 100*tolerance) || (abs(xval-xebs[j+1]) < 100*tolerance && j == nelement) -# ielement = j -# break -# end -# end -# if ielement < 1 -# error("xval=$xval is not within the coordinate $(xcoord.name)") -# end -# return ielement -# end -# end function get_ielement(xval::Float64, xcoord::FiniteElementCoordinate) # search by bisection