Nurbs (#11)

* Add possibility for NURBS Coords

* Test nurbs coords

* fix error from rebase

Author: lijas

src/IgAShell.jl

@@ -9,6 +9,15 @@ using SparseArrays
 using StaticArrays
 using TimerOutputs
 
+"""
+    NurbsCoords and NurbsOrBSplineCoords
+
+Used for dispatching on Nurbs (coords and weights) or Bsplines (only coords), depending on what is used.
+"""
+const NurbsCoords{dim_s,T} = Tuple{Vector{Vec{dim_s,T}}, Vector{T}}
+const NurbsOrBSplineCoords{dim_s,T} = Union{ Vector{Vec{dim_s,T}} , NurbsCoords{dim_s,T}}
+
+
 """
 EdgeInterfaceIndex
     Specifiying interface on cell.

src/igashell_data.jl

@@ -330,7 +330,7 @@ end
 """
 
 """
-struct IGAShellData{dim_p,dim_s,T,LM<:Five.AbstractMaterial,IM<:Five.AbstractCohesiveMaterial}
+struct IGAShellData{dim_p,dim_s,T,LM<:Five.AbstractMaterial,IM<:Five.AbstractCohesiveMaterial,CT}
     layer_materials::Vector{LM}
     interface_material::IM
     orders::NTuple{dim_s,Int}
@@ -351,6 +351,7 @@ struct IGAShellData{dim_p,dim_s,T,LM<:Five.AbstractMaterial,IM<:Five.AbstractCoh
     nqp_ooplane_per_layer::Int
     nqp_interface_order::Int                  
     add_czcells_vtk::Bool
+    coordstype::Type{CT}
 end
 
 function IGAShellData(;
@@ -362,6 +363,7 @@ function IGAShellData(;
     initial_cellstates::Vector{CELLSTATE},
     nqp_inplane_order::Int,
     nqp_ooplane_per_layer::Int,
+    nurbscoords::Bool                       = false,
     width::T                                = 1.0, #Only used in 2d,
     nlayers::Int                            = length(layer_materials),
     zcoords::Vector{T}                      = collect(-thickness/2:(thickness/nlayers):thickness/2),
@@ -384,14 +386,16 @@ function IGAShellData(;
         end
     end
 
+    coordstype = nurbscoords ? NurbsCoords{dim_s,T} : Vec{dim_s, T}
+
     #-----
     ninterfaces = nlayers-1
     ncells = length(initial_cellstates)
     @assert(size(initial_interface_damages)[1] == ninterfaces)
     @assert(size(initial_interface_damages)[2] == ncells)
     dim_s == 3 && @assert(width == 1.0)
 
-    return IGAShellData(layer_materials, interface_material, orders, knot_vectors, thickness, width, nlayers, zcoords, initial_cellstates, initial_interface_damages, adaptable, limit_stress_criterion, limit_damage_criterion, search_radius, collect(locked_elements), small_deformations_theory, nqp_inplane_order, nqp_ooplane_per_layer, nqp_interface_order, add_czcells_vtk)
+    return IGAShellData(layer_materials, interface_material, orders, knot_vectors, thickness, width, nlayers, zcoords, initial_cellstates, initial_interface_damages, adaptable, limit_stress_criterion, limit_damage_criterion, search_radius, collect(locked_elements), small_deformations_theory, nqp_inplane_order, nqp_ooplane_per_layer, nqp_interface_order, add_czcells_vtk, coordstype)
 
 end
 

src/igashell_main.jl

@@ -413,8 +413,10 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
 
     assembler = start_assemble(state.system_arrays.Kⁱ, state.system_arrays.fⁱ, fillzero=false)  
 
-    X = igashell.cache.X
-    Xᵇ = igashell.cache.Xᵇ
+    
+    nnodes = Ferrite.nnodes_per_cell(igashell)
+    X = zeros(igashell.layerdata.coordstype, nnodes) #X = igashell.cache.X
+    Xᵇ = similar(X) # Xᵇ = igashell.cache.Xᵇ
     ue = igashell.cache.ue
     fe = igashell.cache.fe
     ife = igashell.cache.ife
@@ -423,28 +425,26 @@ function _assemble_stiffnessmatrix_and_forcevector!( dh::Ferrite.AbstractDofHand
     ue = igashell.cache.ue
     celldofs = igashell.cache.celldofs
 
-    nqp_oop_per_layer = getnquadpoints_ooplane_per_layer(igashell)
-
     Δt = state.Δt
 
     @timeit "Shell loop" for (ic, cellid) in enumerate(igashell.cellset)
         cellstate = getcellstate(adapdata(igashell), ic)
         cv = build_cellvalue!(igashell, cellstate)
         active_layer_dofs = build_active_layer_dofs(igashell, cellstate)
 
-        Ce = get_extraction_operator(intdata(igashell), ic)
-        IGA.set_bezier_operator!(cv, Ce)
-
         ndofs = Ferrite.ndofs_per_cell(dh, cellid)
         resize!(celldofs, ndofs)
         resize!(ue, ndofs)
 
         Ferrite.cellcoords!(X, dh, cellid)
         Ferrite.celldofs!(celldofs, dh, cellid)
+        
+        Ce = get_extraction_operator(intdata(igashell), ic)
 
         disassemble!(ue, state.d, celldofs)
 
         IGA.compute_bezier_points!(Xᵇ, Ce, X)
+        IGA.set_bezier_operator!(cv, Ce, X)
         @timeit "reinit1" reinit!(cv, Xᵇ)
 
         materialstates = state.partstates[ic].materialstates
@@ -626,9 +626,10 @@ function Five.post_part!(dh, igashell::IGAShell{dim_p,dim_s,T}, states) where {d
         ue = states.d[_celldofs]
 
         nnodes = Ferrite.nnodes_per_cell(igashell)
-        X = zeros(Vec{dim_s,T}, nnodes)
+        X = zeros(igashell.layerdata.coordstype, nnodes)
+        Xᵇ = similar(X)
         Ferrite.cellcoords!(X, dh, cellid)
-        Xᵇ = IGA.compute_bezier_points(Ce, X)
+        IGA.compute_bezier_points!(Xᵇ, Ce, X)
 
         if is_lumped(cellstate)
             _post_lumped(igashell, Xᵇ, X, ue, Ce, cellstate, ic, cellid)
@@ -646,7 +647,7 @@ function _post_layered(igashell, Xᵇ, X, ue, Ce, cellstate, ic::Int, cellid::In
 
     #Shape values for evaluating stresses at center of cell
     cv_mid_interface = igashell.integration_data.cell_value_mid_interfaces
-    set_bezier_operator!(cv_mid_interface, Ce)
+    set_bezier_operator!(cv_mid_interface, Ce, X)
 
     #oop_values = _build_oop_basisvalue!(igashell, cellstate)
     #set_oop_basefunctions!(cv_mid_interface, oop_values)
@@ -678,18 +679,20 @@ function _post_lumped(igashell, Xᵇ, X, ue, Ce, cellstate, ic::Int, cellid::Int
 
     #Build basis_values for cell
     cv = build_cellvalue!(igashell, cellstate)
-    IGA.set_bezier_operator!(cv, Ce)
+    IGA.set_bezier_operator!(cv, Ce, X)
     reinit!(cv, Xᵇ)
 
     #Build basis_values for stress_recovory
     cv_sr = intdata(igashell).cell_values_sr
     oop_values = _build_oop_basisvalue!(igashell, cellstate)
     set_oop_basefunctions!(cv_sr, oop_values)
-    IGA.set_bezier_operator!(cv_sr, Ce)
+    IGA.set_bezier_operator!(cv_sr, Ce, X)
     reinit!(cv_sr, Xᵇ)
 
     recover_cell_stresses(srdata(igashell), σ_states, celldata, cv_sr, cv)
 
+
+    IGA.compute_bezier_points!(Xᵇ, Ce, X)
 end
 
 function _get_layer_forcevector_and_stiffnessmatrix!(

src/igashell_output.jl

@@ -307,6 +307,5 @@ Five.outputname(o::IGAShellConfigStateOutput) = "ElementState"
     dir::Int = 1
 end
 
-
 Five.outputname(o::IGAShellMaterialStateOutput) = string(o.field)
 

src/igashell_utils.jl

@@ -205,4 +205,39 @@ end
 
 #
 _to3d(s::SymmetricTensor{2,2,T}) where T = SymmetricTensor{2,3,T,6}((s[1,1], T(0.0), s[1,2], T(0.0), T(0.0), s[2,2]))
-_to3d(s::SymmetricTensor{2,3}) = s
+_to3d(s::SymmetricTensor{2,3}) = s
+
+"""
+Dispatch on Nurbscoords 
+"""
+
+function Ferrite.cellcoords!((x, w)::NurbsCoords, dh::Ferrite.AbstractDofHandler, i::Int)
+    Ferrite.cellcoords!(x, dh, i)
+    getweights!(w, dh.grid, i)
+end
+
+function IGA.compute_bezier_points!((xb, wb)::NurbsCoords{dim_s,T2}, Ce::BezierExtractionOperator{T}, (x,w)::NurbsCoords{dim_s, T2}; dim::Int=1) where {dim_s,T,T2}
+    IGA.compute_bezier_points!(xb, Ce, w.*x, dim=dim)
+    IGA.compute_bezier_points!(wb, Ce, w)
+    xb .*= inv.(wb)
+end
+
+function IGA.set_bezier_operator!(cv::IGAShellValues, Ce::BezierExtractionOperator{T}, (x,w)::NurbsCoords{dim_s,T}) where {dim_s,T}
+    IGA.set_bezier_operator!(cv, w.*Ce)
+end
+
+function IGA.set_bezier_operator!(cv::IGAShellValues, Ce::BezierExtractionOperator{T}, x::Vector{Vec{dim_s,T}}) where {dim_s,T}
+    IGA.set_bezier_operator!(cv, Ce)
+end
+
+function Base.similar((x,w)::NurbsCoords{dim_s, T}) where {dim_s,T}
+    return similar(x), similar(w)
+end
+
+function Base.zeros(::Type{NurbsCoords{dim_s,T}}, n::Int...) where {dim_s,T}
+    return (zeros(Vec{dim_s,T},n), zeros(T,n))
+end
+
+function myzeros(::Type{NurbsCoords{dim_s,T}}, n::Int) where {dim_s,T}
+    return (zeros(Vec{dim_s,T},n), zeros(T,n))
+end

src/igashell_values.jl

@@ -65,6 +65,7 @@ function function_second_derivative(fe_v::BasisValues{dim_p}, q_point::Int, u::A
     return d²udξ²
 end
 
+Base.copy(cv::BasisValues) = return BasisValues(copy(cv.N), copy(cv.dNdξ), copy(cv.d²Ndξ²))
 
 """
 
@@ -243,22 +244,22 @@ function set_oop_basefunctions!(cv::IGAShellValues{dim_s,dim_p,T}, ooplane_basis
 
 end
 
-function _inplane_nurbs_bezier_extraction(cv::IGAShellValues{dim_s,dim_p,T}, C::IGA.BezierExtractionOperator{T}) where {dim_s,dim_p,T}
-    dBdξ   = cv.inplane_values_bezier.dNdξ
-    B      = cv.inplane_values_bezier.N
+function _inplane_nurbs_bezier_extraction(cv::IGAShellValues{dim_s,dim_p,T}, nurbs::BasisValues, bezier::BasisValues, C::IGA.BezierExtractionOperator{T}) where {dim_s,dim_p,T}
+    dBdξ   = bezier.dNdξ
+    B      = bezier.N
 
     for iq in 1:getnquadpoints_inplane(cv)
         for ib in 1:getnbasefunctions_inplane(cv)
 
-            cv.inplane_values_nurbs.N[ib, iq] = zero(eltype(cv.inplane_values_nurbs.N))
-            cv.inplane_values_nurbs.dNdξ[ib, iq] = zero(eltype(cv.inplane_values_nurbs.dNdξ))
+            nurbs.N[ib, iq] = zero(eltype(nurbs.N))
+            nurbs.dNdξ[ib, iq] = zero(eltype(nurbs.dNdξ))
 
             C_ib = C[ib]
 
             for (i, nz_ind) in enumerate(C_ib.nzind)                
                 val = C_ib.nzval[i]
-                cv.inplane_values_nurbs.N[ib, iq]    += val*   B[nz_ind, iq]
-                cv.inplane_values_nurbs.dNdξ[ib, iq] += val*dBdξ[nz_ind, iq]
+                nurbs.N[ib, iq]    += val*   B[nz_ind, iq]
+                nurbs.dNdξ[ib, iq] += val*dBdξ[nz_ind, iq]
             end
         end
     end
@@ -345,7 +346,7 @@ function _reinit_layer!(cv::IGAShellValues{dim_s,dim_p,T}, qp_indx) where {dim_s
     G[dim_s] = D
 
     detJ = dim_s == 3 ? norm((cross(G[1], G[2]))) : norm(cross(G[1]))
-    cv.detJdA[qp] = detJ*cv.qr.weights[qp]
+    cv.detJdA[qp] = detJ*cv.qr.weights[qp] # TODO: Should be detJ*cv.iqr.weights[iqp] ? 
     cv.detJdV[qp] = detJ*cv.iqr.weights[iqp]*cv.oqr.weights[oqp]*0.5*cv.thickness
 
     Gⁱʲ = inv(SymmetricTensor{2,dim_s,T}((i,j)-> G[i]⋅G[j]))
@@ -420,14 +421,19 @@ function _reinit_midsurface!(cv::IGAShellValues{dim_s,dim_p,T}, iqp::Int, coords
 
 end
 
-function Ferrite.reinit!(cv::IGAShellValues, coords::Vector{Vec{dim_s,T}}) where {dim_s,T}
+function Ferrite.reinit!(cv::IGAShellValues, coords::NurbsOrBSplineCoords{dim_s,T}) where {dim_s,T}
 
     qp = 0
     for iqp in 1:getnquadpoints_inplane(cv)
         _reinit_midsurface!(cv, iqp, coords)
     end
 
-    _inplane_nurbs_bezier_extraction(cv, cv.current_bezier_operator[])
+    if coords isa NurbsCoords
+        _inplane_nurbs_preprocess(cv, coords[2])
+        _inplane_nurbs_bezier_extraction(cv, cv.inplane_values_nurbs, copy(cv.inplane_values_nurbs), cv.current_bezier_operator[])
+    else
+        _inplane_nurbs_bezier_extraction(cv, cv.inplane_values_nurbs, cv.inplane_values_bezier, cv.current_bezier_operator[])
+    end
 
     for oqp in 1:getnquadpoints_ooplane(cv)
         for iqp in 1:getnquadpoints_inplane(cv)
@@ -600,4 +606,123 @@ end
 
 function basis_value(cv::IGAShellValues{dim_s,dim_p,T}, qp::Int, i::Int) where {dim_s,dim_p,T}
     return cv.N[i,qp]
+end
+
+function Ferrite.shape_value(cv::IGAShellValues{dim_s,dim_p,T}, qp::Int, (xᴮ, wᴮ)::NurbsCoords{dim_s,T}) where {dim_s,dim_p,T}
+    val = zero(Vec{dim_s,T})
+    nbasefuncs = getnbasefunctions(cv.inplane_values_bezier)
+    @assert(length(xᴮ) == nbasefuncs)
+
+    W = 0.0
+    @inbounds for i in 1:nbasefuncs
+        W += cv.inplane_values_bezier.N[i,qp] * wᴮ[i]
+    end
+
+    @inbounds for i in 1:nbasefuncs
+        val += xᴮ[i] * wᴮ[i] *inv(W) * cv.inplane_values_bezier.N[i,qp]
+    end
+    return val
+end
+
+function shape_parent_derivative(cv::IGAShellValues{dim_s,dim_p,T}, qp::Int, (xᴮ, wᴮ)::NurbsCoords{dim_s,T}, Θ::Int) where {dim_s,dim_p,T}
+
+    grad = zero(Vec{dim_s,T})
+    nbasefuncs = getnbasefunctions(cv.inplane_values_bezier)
+    @assert(length(xᴮ) == nbasefuncs)
+
+    W = 0.0
+    dWdξ = 0.0
+    @inbounds for i in 1:nbasefuncs
+        N = cv.inplane_values_bezier.N[i,qp]
+        dNdξ = cv.inplane_values_bezier.dNdξ[i,qp][Θ]
+
+        W += N * wᴮ[i]
+        dWdξ +=  dNdξ * wᴮ[i]
+    end
+
+    @inbounds for i in 1:nbasefuncs
+        N = cv.inplane_values_bezier.N[i,qp]
+        dNdξ = cv.inplane_values_bezier.dNdξ[i,qp][Θ]
+        
+        grad += xᴮ[i] * wᴮ[i] * inv(W^2) * (dNdξ * W - dWdξ * N)
+    end
+    return grad
+end
+
+function shape_parent_second_derivative(cv::IGAShellValues{dim_s,dim_p,T}, qp::Int, (xᴮ, wᴮ)::NurbsCoords{dim_s,T}, θ::Tuple{Int,Int}) where {dim_s,dim_p,T}
+    hess = zero(Vec{dim_s,T})
+    nbasefuncs = getnbasefunctions(cv.inplane_values_bezier)
+    @assert(length(xᴮ) == nbasefuncs)
+
+    W = 0.0
+    dWdξ₁ = 0.0
+    dWdξ₂ = 0.0
+    d²Wdξ² = 0.0
+    @inbounds for i in 1:nbasefuncs
+        N = cv.inplane_values_bezier.N[i,qp]
+        dNdξ₁ = cv.inplane_values_bezier.dNdξ[i,qp][θ[1]]
+        dNdξ₂ = cv.inplane_values_bezier.dNdξ[i,qp][θ[2]]
+        d²Ndξ² = cv.inplane_values_bezier.d²Ndξ²[i,qp][θ[1],θ[2]]
+
+        W    += N * wᴮ[i]
+        dWdξ₁ +=  dNdξ₁ * wᴮ[i]
+        dWdξ₂ +=  dNdξ₂ * wᴮ[i]
+        d²Wdξ² +=  d²Ndξ² * wᴮ[i]
+    end
+
+    @inbounds for i in 1:nbasefuncs
+        N = cv.inplane_values_bezier.N[i,qp]
+        dNdξ₁ = cv.inplane_values_bezier.dNdξ[i,qp][θ[1]]
+        dNdξ₂ = cv.inplane_values_bezier.dNdξ[i,qp][θ[2]]
+        d²Ndξ² = cv.inplane_values_bezier.d²Ndξ²[i,qp][θ[1],θ[2]]
+        
+        hess += xᴮ[i] * wᴮ[i] * (
+             -2*inv(W^3) * dWdξ₂ * (W*dNdξ₁ - dWdξ₁*N) + 
+             inv(W^2)*((dWdξ₂*dNdξ₁ + W*d²Ndξ²) - (d²Wdξ²*N + dWdξ₁*dNdξ₂))
+            )
+    end
+    return hess
+end
+
+function Ferrite.spatial_coordinate(cv::IGAShellValues{dim_s,dim_p,T}, qp::Int, x::NurbsCoords{dim_s,T}) where {dim_s,dim_p,T}
+    i2s = CartesianIndices((getnquadpoints_inplane(cv), getnquadpoints_ooplane(cv)))
+    iqp, oqp = Tuple(i2s[qp])
+    D = cv.Eₐ[iqp][dim_s]
+
+    Xᴹ = shape_value(cv, iqp, x)
+    return Xᴹ + cv.oqr.points[oqp][1]*D
+end
+
+function spatial_midsurface_coordinate(cv::IGAShellValues{dim_s,dim_p,T}, iqp::Int, x::NurbsCoords{dim_s,T}) where {dim_s,dim_p,T}
+    @assert(iqp <= getnquadpoints_inplane(cv))
+    Xᴹ = shape_value(cv, iqp, x)
+    return Xᴹ
+end
+
+function _inplane_nurbs_preprocess(cv::IGAShellValues{dim_s,dim_p,T}, wᴮ::Vector{T}) where {dim_s,dim_p,T}
+
+    nbasefuncs = getnbasefunctions(cv.inplane_values_bezier)
+    @assert(length(wᴮ) == nbasefuncs)
+
+    @inbounds for i in 1:getnquadpoints_inplane(cv)
+
+        W = zero(T)
+        dWdξ = zero(Vec{dim_p,T})
+        @inbounds for j in 1:nbasefuncs
+            dNdξ = cv.inplane_values_bezier.dNdξ[j,i]
+            N = cv.inplane_values_bezier.N[j, i]
+
+            W    += wᴮ[j] * N
+            dWdξ += wᴮ[j] * dNdξ
+        end
+
+        @inbounds for j in 1:nbasefuncs
+            dNdξ = cv.inplane_values_bezier.dNdξ[j,i]
+            N = cv.inplane_values_bezier.N[j, i]
+
+            cv.inplane_values_nurbs.dNdξ[j, i] = inv(W)*dNdξ - inv(W^2) * N * dWdξ
+            cv.inplane_values_nurbs.N[j, i] = N/W
+        end
+    end
+
 end