Bisection for get_ielement(), checks on reference grids

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)

src/Interpolation.jl

@@ -1,4 +1,4 @@
-export value_in_coordinate_domain, 
+export value_in_coordinate_domain,
     get_ielement,
     interpolate_1D
 
@@ -8,14 +8,14 @@ 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 
+    else
         in_domain = false
     end
     return in_domain
@@ -27,25 +27,48 @@ assembled grid of the coordinate xcoord.
 """
 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