Merge pull request #73 from sandialabs/safeakima

Safe Akima Interpolation and Extrapolation Fixes

docs/src/literate/C_customizablePreprocessing.jl

@@ -22,7 +22,7 @@ import PyPlot
 import OWENSOpenFASTWrappers
 
 
-path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #splitdir(@__FILE__)[1]
+path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #bld_height_numadsplitdir(@__FILE__)[1]
 
 Inp = OWENS.MasterInput("$runpath/sampleOWENS.yml")
 
@@ -408,11 +408,11 @@ if !AD15On
     delta_xs = shapeX[2:end] - shapeX[1:end-1]
     delta_zs = shapeZ[2:end] - shapeZ[1:end-1]
     delta3D = atan.(delta_xs./delta_zs)
-    delta3D_spl = FLOWMath.akima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
+    delta3D_spl = OWENS.safeakima(shapeZ[1:end-1]./maximum(shapeZ[1:end-1]), delta3D,LinRange(0,1,length(numadIn_bld.span)-1))
 
     bld_height_numad = cumsum(diff(numadIn_bld.span).*(1.0.-abs.(sin.(delta3D_spl)))) # now convert the numad span to a height
 
-    chord = FLOWMath.akima(bld_height_numad./maximum(bld_height_numad), numadIn_bld.chord,LinRange(0,1,Nslices)) # now we can use it to access the numad data 
+    chord = OWENS.safeakima(bld_height_numad./maximum(bld_height_numad), numadIn_bld.chord,LinRange(minimum(bld_height_numad),1,Nslices)) # now we can use it to access the numad data 
     airfoils = fill("nothing",Nslices)
 
     for (iheight_numad,height_numad) in enumerate(bld_height_numad./maximum(bld_height_numad)) # Discretely assign the airfoils
@@ -453,7 +453,7 @@ if AD15On
 
     NumADBldNds = NumADStrutNds = 10 
 
-    bldchord_spl = FLOWMath.akima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
+    bldchord_spl = OWENS.safeakima(numadIn_bld.span./maximum(numadIn_bld.span), numadIn_bld.chord,LinRange(0,1,NumADBldNds))
 
     airfoil_filenames = fill("nothing",NumADBldNds) # Discretely assign the airfoils #TODO: separate out struts
     for (ispan_numad,span_numad) in enumerate(numadIn_bld.span./maximum(numadIn_bld.span))
@@ -470,7 +470,7 @@ if AD15On
         blade_Nnodes = [NumADBldNds for i=1:Nbld]
     else
         blade_filenames = [[blade_filename for i=1:Nbld];[lower_strut_filename for i=1:Nbld];[upper_strut_filename for i=1:Nbld]]
-        strutchord_spl = FLOWMath.akima(numadIn_strut.span./maximum(numadIn_strut.span), numadIn_strut.chord,LinRange(0,1,NumADStrutNds))
+        strutchord_spl = OWENS.safeakima(numadIn_strut.span./maximum(numadIn_strut.span), numadIn_strut.chord,LinRange(0,1,NumADStrutNds))
         blade_chords = [[bldchord_spl for i=1:Nbld];[strutchord_spl for i=1:Nbld];[strutchord_spl for i=1:Nbld]]
         blade_Nnodes = [[NumADBldNds for i=1:Nbld];[NumADStrutNds for i=1:Nbld];[NumADStrutNds for i=1:Nbld]]
     end
@@ -503,7 +503,7 @@ if AD15On
         ymesh = mymesh.y[strt_idx:end_idx]
         ADshapeX = sqrt.(xmesh.^2 .+ ymesh.^2)
         ADshapeX .-= ADshapeX[1] #get it starting at zero #TODO: make robust for blades that don't start at 0
-        ADshapeXspl = FLOWMath.akima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
+        ADshapeXspl = OWENS.safeakima(LinRange(0,H,length(ADshapeX)),ADshapeX,ADshapeZ)
 
         if iADBody<=Nbld #Note that the blades can be curved and are assumed to be oriented vertically
             BlSpn=ADshapeZ

examples/Optimization/optimizationExample.jl

@@ -90,7 +90,7 @@ function messyoptfun!(constraints,Vars)
     x_control_pts_grid = [0.0,0.25,0.5,0.75,1.0]
 
     x_grid = windio[:components][:blade][:internal_structure_2d_fem][:reference_axis][:x][:grid]
-    x_values = FLOWMath.akima(x_control_pts_grid,x_control_pts,x_grid)
+    x_values = safeakima(x_control_pts_grid,x_control_pts,x_grid)
 
     windio[:components][:blade][:internal_structure_2d_fem][:reference_axis][:x][:values] = x_values 
 

examples/SNL34m/34mSetup.jl

@@ -17,7 +17,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 # Analytical
 # z_shape = [0.0, 0.4027099927689326, 0.8054199855378652, 1.2081299783067978, 1.6108399710757304, 2.013549963844663, 2.4162599566135956, 2.818969949382528, 3.221679942151461, 3.624389934920394, 4.027099927689326, 4.429809920458259, 4.832519913227191, 5.235229905996124, 5.637939898765056, 6.221683770581164, 6.821719178571302, 7.437583162221221, 8.068800548394838, 8.714884317956601, 9.375335981533686, 10.049645964128134, 10.737293998282153, 11.437749525493246, 12.305529745531164, 13.201631116890074, 14.122956691386433, 15.066322345375763, 16.028467784161606, 17.00606780965343, 17.995743812332336, 18.994075447807884, 19.997612457611687, 21.002886593374797, 22.006423603178604, 23.004755238654155, 23.99443124133306, 24.97203126682489, 25.934176705610724, 26.877542359600053, 27.79886793409642, 28.694969305455324, 29.562749525493246, 30.263205052704336, 30.950853086858352, 31.62516306945281, 32.285614733029895, 32.93169850259165, 33.56291588876527, 34.1787798724152, 34.77881528040533, 35.362559152221436, 35.821422444858186, 36.280285737494935, 36.739149030131685, 37.198012322768435, 37.656875615405184, 38.11573890804193, 38.57460220067868, 39.03346549331543, 39.49232878595218, 39.951192078588925, 40.410055371225674, 40.868918663862424, 41.32778195649917, 41.78664524913592]
@@ -102,8 +102,8 @@ end
 # PyPlot.figure()
 # PyPlot.plot(bld_OWENSPreCompinput[1].ynode,bld_OWENSPreCompinput[1].ynode)
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_OWENSPreComp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_OWENSPreComp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_OWENSPreComp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_OWENSPreComp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 
@@ -363,12 +363,12 @@ function runowens(model,feamodel,mymesh,myel,aeroForcesDMS,deformTurb;steady=tru
         # samplepts = numadIn_bld.span./maximum(numadIn_bld.span) #normalize #TODO: this is spanwise, while everything else is vertical-wise
         for its = 1:N_ts
             #TODO: there are strain values at each quad point, should be better than just choosing one
-            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#FLOWMath.akima(x,epsilon_x_hist[1,start:stop,its],samplepts)
-            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
-            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#FLOWMath.akima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
-            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_x_hist[1,start:stop,its],samplepts)
-            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_y_hist[1,start:stop,its],samplepts)
-            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#FLOWMath.akima(x,kappa_z_hist[1,start:stop,its],samplepts)
+            eps_x[ibld,its,:] = epsilon_x_hist[1,start:stop,its]#safeakima(x,epsilon_x_hist[1,start:stop,its],samplepts)
+            eps_z[ibld,its,:] = meanepsilon_z_hist[1,start:stop,its]#safeakima(x,meanepsilon_z_hist[1,start:stop,its],samplepts)
+            eps_y[ibld,its,:] = meanepsilon_y_hist[1,start:stop,its]#safeakima(x,meanepsilon_y_hist[1,start:stop,its],samplepts)
+            kappa_x[ibld,its,:] = kappa_x_hist[1,start:stop,its]#safeakima(x,kappa_x_hist[1,start:stop,its],samplepts)
+            kappa_y[ibld,its,:] = kappa_y_hist[1,start:stop,its]#safeakima(x,kappa_y_hist[1,start:stop,its],samplepts)
+            kappa_z[ibld,its,:] = kappa_z_hist[1,start:stop,its]#safeakima(x,kappa_z_hist[1,start:stop,its],samplepts)
         end
     end
 

examples/SNL34m/Example34mHVAWTConfiguration.jl

@@ -81,9 +81,9 @@ SNL34m_5_3_Torque = DelimitedFiles.readdlm("$(path)/data/SAND-91-2228_Data/5.3_T
 
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 offsetTime = 20.0 # seconds
 tocp = [0.0;new_t.+offsetTime; 1e6]
@@ -99,7 +99,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -234,8 +234,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 
 flatwise_stress1grav = (kappa_y_grav[1,end-1,2:end].* thickness .+ 0*eps_x_grav[1,end-1,2:end]) .* Ealuminum
 flatwise_stress2grav = (kappa_y_grav[2,end-1,1:end-1].* thickness .+ 0*eps_x_grav[2,end-1,1:end-1]) .* Ealuminum

examples/SNL34m/Fig4.1Gravity_Fig4.4_Steady40RPMCentrifugalGX.jl

@@ -57,7 +57,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -130,8 +130,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 PyPlot.figure()

examples/SNL34m/Fig5_11_EmergStop_torque_flatwise.jl

@@ -98,12 +98,12 @@ SNL34m_5_11_Vinf = SNL34m_5_11_Vinf[1:2:end,:]
 SNL34m_5_11_RPM[SNL34m_5_11_RPM[:,2].<0.0,2] .*= 0
 
 new_t = LinRange(SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[end,1],100).-SNL34m_5_11_RPM[1,1].+1e-6
-new_RPM = FLOWMath.akima(SNL34m_5_11_RPM[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_11_RPM[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_RPM[:,2],new_t)
 
-new_Torque = FLOWMath.akima(SNL34m_5_11_Torque[:,1],SNL34m_5_11_Torque[:,2],new_t)
+new_Torque = safeakima(SNL34m_5_11_Torque[:,1],SNL34m_5_11_Torque[:,2],new_t)
 
 t_Vinf = LinRange(SNL34m_5_11_Vinf[1,1],SNL34m_5_11_Vinf[end,1],60).-SNL34m_5_11_RPM[1,1].+1e-6
-Vinf_spec = FLOWMath.akima(SNL34m_5_11_Vinf[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_Vinf[:,2],t_Vinf)
+Vinf_spec = safeakima(SNL34m_5_11_Vinf[:,1].-SNL34m_5_11_RPM[1,1],SNL34m_5_11_Vinf[:,2],t_Vinf)
 
 offsetTime = 20.0
 t_Vinf = [0;t_Vinf.+offsetTime;t_Vinf[end]*10]
@@ -127,7 +127,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -200,8 +200,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]
 
 PyPlot.figure()

examples/SNL34m/Fig5_4_Fig5_3_normaloperation_torque_flatwise_generatorControl.jl

@@ -78,11 +78,11 @@ PyPlot.xlim([0,100.0])
 # PyPlot.savefig("$(path)/../figs/NormalOperation_Vinf.pdf",transparent = true)
 
 new_t = LinRange(SNL34m_5_3_RPM[1,1],SNL34m_5_3_RPM[end,1],100)
-new_RPM = FLOWMath.akima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
+new_RPM = safeakima(SNL34m_5_3_RPM[:,1],SNL34m_5_3_RPM[:,2],new_t)
 
-new_Torque = FLOWMath.akima(SNL34m_5_3_Torque[:,1],SNL34m_5_3_Torque[:,2],new_t)
+new_Torque = safeakima(SNL34m_5_3_Torque[:,1],SNL34m_5_3_Torque[:,2],new_t)
 
-Vinf_spec = FLOWMath.akima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
+Vinf_spec = safeakima(SNL34m_5_3_Vinf[:,1],SNL34m_5_3_Vinf[:,2],new_t)
 
 t_Vinf = [0;new_t;1e6]
 Vinf_spec = [Vinf_spec[1];Vinf_spec;Vinf_spec[end]]
@@ -121,7 +121,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = FLOWMath.akima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -194,8 +194,8 @@ bld1start = Int(mymesh.structuralNodeNumbers[1,1])
 bld1end = Int(mymesh.structuralNodeNumbers[1,end])
 spanpos = [0.0;cumsum(sqrt.(diff(mymesh.x[bld1start:bld1end]).^2 .+ diff(mymesh.z[bld1start:bld1end]).^2))]
 spanposmid = cumsum(diff(spanpos))
-thickness = FLOWMath.akima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
-thickness_lag = FLOWMath.akima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
+thickness = safeakima(numadIn_bld.span,thickness_precomp_flap,spanposmid)
+thickness_lag = safeakima(numadIn_bld.span,thickness_precomp_lag,spanposmid)
 # thickness = thicknessGX[1:end-1]