remove conncectivy matrix, and use cellnodes!(dh) instead

Author: lijas

examples/enf_2d.jl

@@ -86,7 +86,6 @@ IgAShell.IGAShellData(;
 igashell = 
 IgAShell.IGAShell(
     cellset = partset1, 
-    connectivity = nurbsmesh.IEN,#reverse(nurbsmesh.IEN, dims=1), 
     data = igashelldata
 ) 
 push!(data.parts, igashell)

src/igashell_adaptivity.jl

@@ -8,6 +8,7 @@ struct IGAShellAdaptivity{T}
     cellstates::Vector{CELLSTATE}
     control_point_states::Vector{CPSTATE}
     propagation_checked::Matrix{Bool}
+    locked_elements::Vector{Int}
     locked_control_points::Vector{Int}
 
     lumped2layered::IGA.BezierExtractionOperator{T}
@@ -32,7 +33,6 @@ function get_controlpoint_state(adap::IGAShellAdaptivity, cpid::Vector{Int})
     return adap.control_point_states[cpid]
 end
 
-
 function set_controlpoint_state!(adap::IGAShellAdaptivity, cpid::Int, state::CPSTATE)
     adap.control_point_states[cpid] = state
 end
@@ -41,7 +41,7 @@ function setcellstate!(adap::IGAShellAdaptivity, cellid::Int, state::CELLSTATE)
     adap.cellstates[cellid] = state
 end
 
-function IGAShellAdaptivity(data::IGAShellData{dim_p,dim_s,T}, cell_connectivity::Matrix{Int}, ncells, nnodes) where {dim_p,dim_s,T}
+function IGAShellAdaptivity(data::IGAShellData{dim_p,dim_s,T}, ncells) where {dim_p,dim_s,T}
 
     order = data.orders[dim_s]
     ninterfaces = data.nlayers-1
@@ -66,49 +66,58 @@ function IGAShellAdaptivity(data::IGAShellData{dim_p,dim_s,T}, cell_connectivity
     Cla2di = IGA.bezier_extraction_to_vector(Cmat_la2di')
     Clu2di = IGA.bezier_extraction_to_vector(Cmat_lu2di')
 
+    propagation_checked = [false for _ in 1:ninterfaces, _ in 1:ncells]
+    
+    node_states = CPSTATE[] #initialized later in _init_cpstates_cellstates
+    locked_control_points = Int[]
+
+    return IGAShellAdaptivity(data.initial_cellstates, node_states, propagation_checked, data.locked_elements, locked_control_points,
+                              Clu2la, Cla2di, Clu2di, 
+                              weakdiscont2discont, strongdiscont2discont,
+                              interface_knots, order)
+end
+
+function _init_cpstates_cellstates!(dh::Ferrite.AbstractDofHandler, igashell)
+
+    nnodes = get_n_controlpoints(igashell.layerdata)
+    cellnodes = igashell.cache.cellnodes
+
     #Some of the cells will be initialized with LUMPED, and some with LYARED/DISCONTINIUOS
     # This means that some cells will be in a mixed mode... Determine those
     # Prioritize the DISCONTINIUOS before LAYERED and LUMPED, ie. if one node is both lumped 
     # and disocontinous, choose it to be dinscontionous
     node_states = fill(LUMPED_CPSTATE, nnodes)
 
-    for cellid in 1:ncells
-        cellstate = data.initial_cellstates[cellid] 
-        for cellnodes in cell_connectivity[:, cellid]
-            for nodeid in cellnodes
-                node_states[nodeid] = combine_states(node_states[nodeid], first(cellstate.cpstates), ninterfaces)
-            end
+    for cellid in 1:getncells(igashell)
+        cellstate = igashell.adaptivity.cellstates[cellid] 
+        Ferrite.cellnodes!(cellnodes, dh, cellid)
+        for nodeid in cellnodes
+            node_states[nodeid] = combine_states(node_states[nodeid], first(cellstate.cpstates), ninterfaces(igashell))
         end
+
     end
 
     #Check if there is any cell that is has MIXED states controlpoints
-    for cellid in 1:ncells
-        nodeids = cell_connectivity[:, cellid]
-        cellnode_states = node_states[nodeids]
+    for cellid in 1:getncells(igashell)
+        Ferrite.cellnodes!(cellnodes, dh, cellid)
+        cellnode_states = node_states[cellnodes]
 
         #if NOT all of the nodes in the cell are equal, the element is mixed
         if !all(first(cellnode_states) .== cellnode_states)
-            data.initial_cellstates[cellid] = CELLSTATE(_MIXED, copy(cellnode_states))
+            igashell.adaptivity.cellstates[cellid] = CELLSTATE(_MIXED, copy(cellnode_states))
         end
 
     end
-    
-    propagation_checked = [false for _ in 1:ninterfaces, _ in 1:ncells]
 
     #Some cells are not allowed to be upgraded (locked).. store the contorl points of these cells
-    locked_control_points = Int[]
-    for cellid in data.locked_elements
-        for cell_nodes in cell_connectivity[:, cellid]
-            for nodeid in cell_nodes
-                push!(locked_control_points, nodeid)
-            end
+    for cellid in igashell.adaptivity.locked_elements
+        Ferrite.cellnodes!(cellnodes, dh, cellid)
+        for nodeid in cellnodes
+            push!(igashell.adaptivity.locked_control_points, nodeid)
         end
     end
 
-    return IGAShellAdaptivity(data.initial_cellstates, node_states, propagation_checked, locked_control_points,
-                              Clu2la, Cla2di, Clu2di, 
-                              weakdiscont2discont, strongdiscont2discont,
-                              interface_knots, order)
+    Ferrite.copy!!(igashell.adaptivity.control_point_states, node_states)
 end
 
 function get_upgrade_operator(adap::IGAShellAdaptivity, from::CPSTATE, to::CPSTATE)

src/igashell_data.jl

@@ -417,4 +417,5 @@ getwidth(data::IGAShellData) = data.width
 
 LIMIT_STRESS_CRITERION(data::IGAShellData) = data.limit_stress_criterion
 LIMIT_DAMAGE_VARIABLE(data::IGAShellData) = data.limit_damage_criterion
-PROPAGATION_SEARCH_RADIUS(data::IGAShellData) = data.search_radius
+PROPAGATION_SEARCH_RADIUS(data::IGAShellData) = data.search_radius
+get_n_controlpoints(data::IGAShellData{dim_p}) where {dim_p} = prod( (i) -> length(data.knot_vectors[i]) - data.orders[i] - 1, 1:dim_p)

src/igashell_main.jl

@@ -1,5 +1,13 @@
 export IGAShellState, IGAShell
 
+struct IGAShellCache
+    cellnodes::Vector{Int}
+    #celldofs::Vector{Int}
+    #cellcoords::Vector{Vec{dim,T}}
+end
+
+IGAShellCache(data::IGAShellData) = IGAShellCache(zeros(Int, Ferrite.nnodes_per_cell(data)))
+
 """
     IGAShell
 
@@ -10,16 +18,17 @@ struct IGAShell{dim_p, dim_s, T,
                 intdata<:IGAShellIntegrationData, 
                 adapdata<:IGAShellAdaptivity, 
                 vtkdata<:IGAShellVTK,
-                srdata<:IGAShellStressRecovory} <: Five.AbstractPart{dim_s}
+                srdata<:IGAShellStressRecovory,
+                cache<:IGAShellCache} <: Five.AbstractPart{dim_s}
 
     layerdata::data
     integration_data::intdata
     adaptivity::adapdata
     vtkdata::vtkdata
     stress_recovory::srdata
+    cache::cache
 
     cellset::Vector{Int}
-    cell_connectivity::Matrix{Int}
 end
 
 #Utility functions
@@ -29,9 +38,6 @@ intdata(igashell::IGAShell) = igashell.integration_data
 vtkdata(igashell::IGAShell) = igashell.vtkdata
 srdata(igashell::IGAShell) = igashell.stress_recovory
 
-cellconectivity!(nodes::Vector{Int}, igashell::IGAShell, cellid::Int) =  nodes .= igashell.cell_connectivity[:, cellid]
-cellconectivity(igashell::IGAShell, cellid::Int) = @view igashell.cell_connectivity[:, cellid]
-
 Ferrite.getnquadpoints(igashell::IGAShell) = getnquadpoints_per_layer(igashell)*nlayers(igashell)
 getnquadpoints_inplane(igashell::IGAShell) = length(getweights(intdata(igashell).iqr))
 getnquadpoints_ooplane(igashell::IGAShell) = length(getweights(intdata(igashell).oqr))
@@ -61,7 +67,6 @@ getcellstate(igashell::IGAShell, i::Int) = igashell.adaptivity.cellstates[i]
 Ferrite.nnodes_per_cell(igashell::IGAShell{dim_p}, cellid::Int=1) where dim_p = prod(igashell.layerdata.orders[1:dim_p].+1)::Int#getnbasefunctions(igashell.cv_inplane) ÷ dim_p
 Ferrite.getdim(igashell::IGAShell{dim_p,dim_s}) where {dim_p,dim_s} = dim_s
 Ferrite.getncells(igashell::IGAShell) = length(igashell.cellset)
-Ferrite.getnnodes(igashell::IGAShell) = maximum(igashell.cell_connectivity)
 
 Five.get_fields(igashell::IGAShell) = [Field(:u, getmidsurface_ip(layerdata(igashell)), ndofs_per_controlpoint(igashell, LUMPED_CPSTATE))]
 
@@ -74,7 +79,7 @@ function get_inplane_qp_range(n_qp_per_layer::Int, n_inpqp::Int, ilayer::Int, ro
     return (1:n_inpqp) .+ offset
 end
 
-function _igashell_input_checks(data::IGAShellData{dim_p, dim_s}, cellset::AbstractVector{Int}, cell_connectivity::Matrix{Int}) where {dim_s,dim_p}
+function _igashell_input_checks(data::IGAShellData{dim_p, dim_s}, cellset::AbstractVector{Int}) where {dim_s,dim_p}
 
     @assert(!any(is_mixed.(data.initial_cellstates)))
     @assert( dim_s == length(data.orders) )
@@ -83,17 +88,15 @@ end
 
 function IGAShell(;
     cellset::AbstractVector{Int}, 
-    connectivity::Matrix{Int},
     data::IGAShellData{dim_p,dim_s,T}
     ) where {dim_p,dim_s,T}
 
-    _igashell_input_checks(data, cellset, connectivity)
+    _igashell_input_checks(data, cellset)
 
     ncells = length(cellset)
-    ncontrol_points = maximum(connectivity)
 
     #Setup adaptivity structure
-    adapdata = IGAShellAdaptivity(data, connectivity, ncells, ncontrol_points)
+    adapdata = IGAShellAdaptivity(data, ncells)
 
     #
     Ce_mat, _ = IGA.compute_bezier_extraction_operators(data.orders[1:dim_p], data.knot_vectors[1:dim_p])
@@ -103,13 +106,15 @@ function IGAShell(;
     vtkdata = IGAShellVTK(data)
     intdata = IGAShellIntegrationData(data, Ce_vec)
     srdata = IGAShellStressRecovory(data)
+    cache = IGAShellCache(data)
 
-    return IGAShell{dim_p,dim_s,T, typeof(data), typeof(intdata), typeof(adapdata), typeof(vtkdata), typeof(srdata)}(data, intdata, adapdata, vtkdata, srdata, cellset, connectivity)
+    return IGAShell{dim_p,dim_s,T, typeof(data), typeof(intdata), typeof(adapdata), typeof(vtkdata), typeof(srdata), typeof(cache)}(data, intdata, adapdata, vtkdata, srdata, cache, cellset)
 
 end
 
 function Five.init_part!(igashell::IGAShell, dh::Ferrite.AbstractDofHandler)
     _init_vtk_grid!(dh, igashell)
+    _init_cpstates_cellstates!(dh, igashell)
 end
 
 struct IGAShellState{MS1<:Five.AbstractMaterialState,MS2<:Five.AbstractMaterialState} <: Five.AbstractPartState

src/igashell_upgrade.jl

@@ -5,7 +5,7 @@ function _commit_part!(dh::Ferrite.AbstractDofHandler,
 
     instructions = FieldDimUpgradeInstruction[]
     upgraded_cells = Int[]
-
+    cellnodes = zeros(Int, Ferrite.nnodes_per_cell(igashell))
     #copy of the current node states
     node_states = deepcopy(adapdata(igashell).control_point_states)
 
@@ -61,7 +61,8 @@ function _commit_part!(dh::Ferrite.AbstractDofHandler,
 
             # Loop through all nodes for this cell and set the state 
             # of the node to the state of the cell, if it is an "upgrade"
-            for (i, nodeid) in enumerate(igashell.cell_connectivity[:, ic])
+            Ferrite.cellnodes!(cellnodes, dh, cellid)
+            for (i, nodeid) in enumerate(cellnodes)
                 node_states[nodeid] = combine_states(node_states[nodeid], get_cpstate(upgrade_to, i), ninterfaces(igashell))
             end
         end
@@ -76,7 +77,7 @@ function _commit_part!(dh::Ferrite.AbstractDofHandler,
         ue = state.d[_celldofs]
         Δue = state.Δd[_celldofs]
 
-        cellnodes = igashell.cell_connectivity[:, ic]
+        Ferrite.cellnodes!(cellnodes, dh, cellid)
         new_cellnode_states = node_states[cellnodes]
         current_cellnode_states = igashell.adaptivity.control_point_states[cellnodes]
 
@@ -86,7 +87,7 @@ function _commit_part!(dh::Ferrite.AbstractDofHandler,
             continue
         end
 
-        instr = construct_upgrade_instruction(igashell, cellid, current_cellnode_states, new_cellnode_states, ue, Δue)
+        instr = construct_upgrade_instruction(igashell, cellid, cellnodes, current_cellnode_states, new_cellnode_states, ue, Δue)
         push!(instructions, instr)
 
         # Check if the element is mixed.
@@ -236,7 +237,7 @@ function stress_interface_index(qplayerid::Int, nqpinplane::Int, dim_p::Int)
     return idx
 end
 
-function construct_upgrade_instruction(igashell::IGAShell{dim_p,dim_s,T}, cellid::Int, current_cellnode_state::AbstractVector{CPSTATE}, cellnode_states::AbstractVector{CPSTATE}, ue::Vector{T}, Δue::Vector{T}) where {dim_p,dim_s,T}
+function construct_upgrade_instruction(igashell::IGAShell{dim_p,dim_s,T}, cellid::Int, cellnodes::Vector{Int}, current_cellnode_state::AbstractVector{CPSTATE}, cellnode_states::AbstractVector{CPSTATE}, ue::Vector{T}, Δue::Vector{T}) where {dim_p,dim_s,T}
 
     local_dof_idxs = Int[]
     current_fielddims = Int[]
@@ -250,7 +251,7 @@ function construct_upgrade_instruction(igashell::IGAShell{dim_p,dim_s,T}, cellid
     # Loop through all node in the cell and create the "instructions" that the 
     # dofhandler needs in order to distribute new dofs.
 
-    for (i, nodeid) in enumerate(igashell.cell_connectivity[:,cellid])
+    for (i, nodeid) in enumerate(cellnodes)
 
         cp_state = current_cellnode_state[i]
         current_fielddim = ndofs_per_controlpoint(igashell, cp_state)

src/igashell_utils.jl

@@ -37,18 +37,19 @@ end
 function initial_upgrade_of_dofhandler(dh::MixedDofHandler, igashell::IGAShell)
 
     instructions = Five.FieldDimUpgradeInstruction[]
+    cellnodes = igashell.cache.cellnodes
 
     for (ic, cellid) in enumerate(igashell.cellset)
 
-        cellnodes = igashell.cell_connectivity[:, ic]
+        Ferrite.cellnodes!(cellnodes, dh, cellid)
         cellnode_states = adapdata(igashell).control_point_states[cellnodes]
 
         initial_cellnode_states = fill(LUMPED_CPSTATE, length(cellnode_states))
 
         if cellnode_states != initial_cellnode_states
             ndofs = ndofs_per_cell(dh, cellid)
 
-            instr = construct_upgrade_instruction(igashell, cellid, initial_cellnode_states, cellnode_states, zeros(Float64, ndofs), zeros(Float64, ndofs))
+            instr = construct_upgrade_instruction(igashell, cellid, cellnodes, initial_cellnode_states, cellnode_states, zeros(Float64, ndofs), zeros(Float64, ndofs))
             push!(instructions, instr)
         end
 
@@ -135,34 +136,6 @@ function active_basefunctions(field_dim::Int, ::Int)
     return 1:field_dim
 end
 
-
-#igashelldofs(dh::Ferrite.AbstractDofHandler, igashell::IGAShell, index::GeometryObject) = celldofs(dh, index[1])[igashelldofs(igashell, index)]
-
-function igashelldofs(igashell::IGAShell{dim_p,dim_s}, index::GeometryObject, components::Vector{Int}=collect(1:dim_s)) where {dim_p,dim_s}
-
-    cellid = index[1]
-    ip = getmidsurface_ip(layerdata(igashell))
-
-    facepoints = geometryobject(ip, index)
-
-    face_dofs = Int[]; currentdof = 1
-    for (i, nodeid) in enumerate(cellconectivity(igashell, cellid))
-        nodestate = get_controlpoint_state(adapdata(igashell), nodeid)
-        nnodes_per_controlpoints = ndofs_per_controlpoint(igashell, nodestate) ÷ dim_s
-        if i in facepoints
-            for basefunc in active_basefunctions(nnodes_per_controlpoints, index)
-                for d in components
-                    push!(face_dofs, currentdof + (basefunc-1)*dim_s + d -1)
-                end
-            end
-        end
-        currentdof += nnodes_per_controlpoints*dim_s
-    end
-
-    return face_dofs
-end
-
-
 function extrapolate(y1,y2,x1,x2,x)
     return y1 + ((x-x1)/(x2-x1))*(y2-y1)
 end#

src/igashell_vtk.jl

@@ -162,7 +162,7 @@ function Five.get_vtk_displacements(dh::Ferrite.AbstractDofHandler, igashell::IG
         ue = state.d[celldofs]
 
         #Transform displayements to fully discontinuous state
-        cellconectivity!(nodes, igashell, ic)
+        Ferrite.cellnodes!(nodes, dh, ic)
         ue_bezier = T[]
         offset = 1
 
@@ -211,14 +211,16 @@ function Five.get_vtk_celldata(igashell::IGAShell{dim_p,dim_s,T}, output::VTKCel
     end
 
     cellstates = zeros(Int, nvtkcells)
+    cellnodes = zeros(Int, Ferrite.nnodes_per_cell(igashell))
+
     for (ic, cellid) in enumerate(igashell.cellset)
         cellstate = getcellstate(adapdata(igashell), ic)
         for ilay in 1:nlayers(igashell)
             for vtkcell in nvtkcells_per_layer(vtkdata(igashell))
-                cellnodes = zeros(Int, Ferrite.nnodes_per_cell(igashell, ic))
-                cellconectivity = cellconectivity!(cellnodes, igashell, ic)
-                tmp_cellstate = get_controlpoint_state(adapdata(igashell), cellconectivity[1])
-                for (i, nodeid) in enumerate(cellconectivity[2:end])
+                
+                Ferrite.cellnodes!(cellnodes, globaldata.dh, ic)
+                tmp_cellstate = get_controlpoint_state(adapdata(igashell), cellnodes[1])
+                for (i, nodeid) in enumerate(cellnodes[2:end])
                     cp_state = get_controlpoint_state(adapdata(igashell), nodeid)
                     tmp_cellstate = tmp_cellstate.state > cp_state.state ? tmp_cellstate : cp_state
                 end

test/igashell_test.jl

@@ -99,7 +99,6 @@ function get_curved_mesh(cellstate; h, b, R)
 
     igashell = IgAShell.IGAShell(
         cellset = collect(1:getncells(grid)), 
-        connectivity = reverse(nurbsmesh.IEN, dims=1), 
         data = igashelldata) 
 
     return grid, igashell
@@ -195,7 +194,6 @@ function get_test_mesh(CELLSTATE::IgAShell.CELLSTATE, nelx, damage)
     igashell = 
     IgAShell.IGAShell(
         cellset = 1:nelx, 
-        connectivity = reverse(nurbsmesh.IEN, dims=1), 
         data = igashelldata
     ) 
     push!(data.parts, igashell)