Merge pull request #89 from byuflowlab/no_webs

Fix bug for when the layup contains no webs - closes #88

src/afmesh.jl

@@ -20,6 +20,9 @@ end
 Base.eltype(::Layer{TF}) where {TF} = TF
 Base.eltype(::Type{Layer{TF}}) where {TF} = TF
 
+Layer{TF}(l::Layer) where {TF} = Layer{TF}(l.material, l.t, l.theta)
+Base.convert(::Type{Layer{TF}}, l::Layer) where {TF} = Layer{TF}(l)
+
 # from https://discourse.julialang.org/t/findnearest-function/4143/4
 function searchsortednearest(a, x)
     idx = searchsortedfirst(a, x)
@@ -848,7 +851,9 @@ function afmesh(xaf, yaf, chord, twist, paxis, xbreak, webloc, segments, webs; d
     nlayer = length(segments[1])
     nodes, elements = combine_halfs(nodesu, elementsu, nodesl, elementsl, nlayer, x_te)
 
-    nodes, elements = addwebs(idx_webu, idx_webl, nx_web, nodes, elements, webs, length(nodesu), nlayer, wns)
+    if nw > 0 # only add webs if there are webs defined
+        nodes, elements = addwebs(idx_webu, idx_webl, nx_web, nodes, elements, webs, length(nodesu), nlayer, wns)
+    end
     # -----------------------------------
 
     # ------ rotate with twist -------

src/section.jl

@@ -27,6 +27,10 @@ end
 Base.eltype(::Material{TF}) where TF = TF
 Base.eltype(::Type{Material{TF}}) where TF = TF
 
+Material{TF}(m::Material) where {TF} = Material{TF}(m.E1, m.E2, m.E3, m.G12, m.G13, m.G23, 
+    m.nu12, m.nu13, m.nu23, m.rho)
+Base.convert(::Type{Material{TF}}, m::Material) where {TF} = Material{TF}(m)
+
 """
     Node(x, y)
 
@@ -44,6 +48,9 @@ end
 Base.eltype(::Node{TF}) where TF = TF
 Base.eltype(::Type{Node{TF}}) where TF = TF
 
+Node{TF}(n::Node) where {TF} = Node{TF}(n.x, n.y)
+Base.convert(::Type{Node{TF}}, n::Node) where {TF} = Node{TF}(n)
+
 """
     MeshElement(nodenum, material, theta)
 
@@ -64,6 +71,9 @@ end
 Base.eltype(::MeshElement{VI, TF}) where {VI, TF} = TF
 Base.eltype(::Type{MeshElement{VI, TF}}) where {VI, TF} = TF
 
+MeshElement{VI,TF}(e::MeshElement) where {VI,TF} = MeshElement{VI,TF}(e.nodenum, e.material, e.theta)
+Base.convert(::Type{MeshElement{VI,TF}}, e::MeshElement) where {VI,TF} = MeshElement{VI,TF}(e)
+
 """
 internal cache so allocations happen only once upfront
 """

test/section.jl

@@ -713,7 +713,6 @@ function sectionwrapper(x)
 
     webs = [web, web]
 
-
     nodes, elements = afmesh(xaf, yaf, chord, twist, paxis, xbreak, webloc, segments, webs)
 
     S, sc, tc = compliance_matrix(nodes, elements)
@@ -722,6 +721,57 @@ function sectionwrapper(x)
     return vcat([S; M]...)
 end
 
+function sectionwrapper_nowebs(x)
+
+    TF = eltype(x)
+
+    xaf = TF[1.00000000, 0.99619582, 0.98515158, 0.96764209, 0.94421447, 0.91510964, 0.88074158, 0.84177999, 0.79894110, 0.75297076, 0.70461763, 0.65461515, 0.60366461, 0.55242353, 0.50149950, 0.45144530, 0.40276150, 0.35589801, 0.31131449, 0.26917194, 0.22927064, 0.19167283, 0.15672257, 0.12469599, 0.09585870, 0.07046974, 0.04874337, 0.03081405, 0.01681379, 0.00687971, 0.00143518, 0.00053606, 0.00006572, 0.00001249, 0.00023032, 0.00079945, 0.00170287, 0.00354717, 0.00592084, 0.01810144, 0.03471169, 0.05589286, 0.08132751, 0.11073805, 0.14391397, 0.18067874, 0.22089879, 0.26433734, 0.31062190, 0.35933893, 0.40999990, 0.46204424, 0.51483073, 0.56767889, 0.61998250, 0.67114514, 0.72054815, 0.76758733, 0.81168064, 0.85227225, 0.88883823, 0.92088961, 0.94797259, 0.96977487, 0.98607009, 0.99640466, 1.00000000]
+    yaf = TF[0.00000000, 0.00017047, 0.00100213, 0.00285474, 0.00556001, 0.00906779, 0.01357364, 0.01916802, 0.02580144, 0.03334313, 0.04158593, 0.05026338, 0.05906756, 0.06766426, 0.07571157, 0.08287416, 0.08882939, 0.09329359, 0.09592864, 0.09626763, 0.09424396, 0.09023579, 0.08451656, 0.07727756, 0.06875796, 0.05918984, 0.04880096, 0.03786904, 0.02676332, 0.01592385, 0.00647946, 0.00370956, 0.00112514, -0.00046881, -0.00191488, -0.00329201, -0.00470585, -0.00688469, -0.00912202, -0.01720842, -0.02488211, -0.03226730, -0.03908459, -0.04503763, -0.04986836, -0.05338180, -0.05551392, -0.05636585, -0.05605816, -0.05472399, -0.05254383, -0.04969990, -0.04637175, -0.04264894, -0.03859653, -0.03433153, -0.02996944, -0.02560890, -0.02134397, -0.01726049, -0.01343567, -0.00993849, -0.00679919, -0.00402321, -0.00180118, -0.00044469, 0.00000000]
+
+    xaf[15] = x[1]
+    yaf[15] = x[2]
+
+    uni = Material{TF}(37.00e9, 9.00e9, 9.00e9, 4.00e9, 4.00e9, 4.00e9, 0.28, 0.28, 0.28, 1.86e3)
+    double = Material{TF}(x[3], 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.83e3)
+    gelcoat = Material{TF}(1e1, 1e1, 1e1, 1.0, 1.0, 1.0, 0.30, 0.30, 0.30, 1.83e3)
+    nexus = Material{TF}(10.30e9, 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.664e3)
+    balsa = Material{TF}(0.01e9, 0.01e9, 0.01e9, 2e5, 2e5, 2e5, 0.30, 0.30, 0.30, 0.128e3)
+    mat = [uni, double, gelcoat, nexus, balsa]
+
+    chord = x[4]
+    twist = x[5]
+    paxis = 0.4750 / chord
+    xbreak = [0.0, 0.0041, 0.1147, 0.5366, 1.0]
+    webloc = []
+
+    idx = [3, 4, 2]
+    t = TF[0.000381, 0.00051, 18*0.00053]
+    theta = TF[0, 0, 20]*pi/180
+    layup1 = Layer.(mat[idx], t, theta)
+    idx = [3, 4, 2]
+    t = TF[0.000381, 0.00051, 33*0.00053]
+    theta = TF[0, 0, 20]*pi/180
+    layup2 = Layer.(mat[idx], t, theta)
+    idx = [3, 4, 2, 1, 5, 1, 2]
+    t = TF[0.000381, 0.00051, x[6], 38*0.00053, 1*0.003125, 37*0.00053, 16*0.00053]
+    theta = TF[0, 0, x[7], 30, 0, 30, 20]*pi/180
+    layup3 = Layer.(mat[idx], t, theta)
+    idx = [3, 4, 2, 5, 2]
+    t = TF[0.000381, 0.00051, 17*0.00053, 0.003125, 16*0.00053]
+    theta = TF[0, 0, 20, 0, 0]*pi/180
+    layup4 = Layer.(mat[idx], t, theta)
+
+    segments = [layup1, layup2, layup3, layup4]
+
+    webs = []
+
+    nodes, elements = afmesh(xaf, yaf, chord, twist, paxis, xbreak, webloc, segments, webs)
+
+    S, sc, tc = compliance_matrix(nodes, elements)
+    M, mc = mass_matrix(nodes, elements)
+
+    return vcat([S; M]...)
+end
 
 @testset "type stability" begin
 
@@ -740,7 +790,7 @@ end
 
 end
 
-@testset "Jacobian" begin
+@testset "Jacobian with Webs" begin
 
     # should use graph coloring b.c. plenty of sparsity, but dense is fine for purpose of this test.
     xvec = [0.5014995, 0.07571157, 10.30e9, 1.9, 0.0*pi/180, 17*0.00053, 20.0, 38*0.00053, 0.0]
@@ -751,3 +801,14 @@ end
     @test maximum(abs.(J1 .- J2)) < 1e-6
 
 end
+
+@testset "Jacobian without Webs" begin
+
+    xvec = [0.5014995, 0.07571157, 10.30e9, 1.9, 0.0*pi/180, 17*0.00053, 20.0, 38*0.00053, 0.0]
+
+    J1 = ForwardDiff.jacobian(sectionwrapper_nowebs, xvec)
+    J2 = FiniteDiff.finite_difference_jacobian(sectionwrapper_nowebs, xvec, Val{:central})
+
+    @test maximum(abs.(J1 .- J2)) < 1e-6
+
+end