Merge pull request #234 from kaipartmann/bc_handling

Introducing a boundary condition handler

Author: kaipartmann

docs/src/private_api_reference.md

@@ -23,6 +23,7 @@ Peridynamics.find_points
 Peridynamics.apply_precracks!
 Peridynamics.apply_precrack!
 Peridynamics.point_sets_intersect
+Peridynamics.displacement_bc!
 Peridynamics.velocity_databc!
 Peridynamics.forcedensity_databc!
 Peridynamics.invreg
@@ -48,6 +49,7 @@ Peridynamics.InteractionSystem
 Peridynamics.PointSetsPreCrack
 Peridynamics.StandardPointParameters
 Peridynamics.SingleDimBC
+Peridynamics.PosSingleDimBC
 Peridynamics.PosDepSingleDimBC
 Peridynamics.CKIPointParameters
 Peridynamics.BACPointParameters

src/Peridynamics.jl

@@ -116,6 +116,7 @@ include("discretization/bond_system_corrections.jl")
 include("discretization/zem_stabilization.jl")
 include("discretization/bond_associated_system.jl")
 include("discretization/interaction_system.jl")
+include("discretization/condition_handler.jl")
 include("discretization/body_chunk.jl")
 
 include("core/job.jl")

src/auxiliary/io.jl

@@ -167,5 +167,5 @@ function export_fields!(vtk, chunk, fields_spec::Dict{Symbol,Vector{Symbol}})
 end
 
 @inline function get_loc_position(chunk::AbstractBodyChunk)
-    return @views chunk.storage.position[:, 1:n_loc_points(chunk)]
+    return @views chunk.storage.position[:, 1:get_n_loc_points(chunk)]
 end

src/core/systems.jl

@@ -16,6 +16,10 @@ end
     return get_n_loc_points(system.chunk_handler)
 end
 
+@inline function get_n_points(system::AbstractSystem)
+    return get_n_points(system.chunk_handler)
+end
+
 @inline function get_point_ids(system::AbstractSystem)
     return get_point_ids(system.chunk_handler)
 end
@@ -48,10 +52,78 @@ end
     return get_loc_view(a, system.chunk_handler)
 end
 
-@inline function get_n_points(system::AbstractSystem)
-    return get_n_points(system.chunk_handler)
+@inline function get_n_dim(::AbstractSystem)
+    return 3 # 3D only for now
+end
+
+@inline function get_n_dof(system::AbstractSystem)
+    return get_n_dim(system) * get_n_points(system)
+end
+
+@inline function get_n_loc_dof(system::AbstractSystem)
+    return get_n_dim(system) * get_n_loc_points(system)
+end
+
+@inline function each_dim(system::AbstractSystem)
+    return 1:get_n_dim(system)
+end
+
+@inline function get_dof(system::AbstractSystem, dim::Int, point_id::Int)
+    return get_dof(get_n_dim(system), dim, point_id)
+end
+
+@inline function get_dof(n_dim::Int, dim::Int, point_id::Int)
+    return (point_id - 1) * n_dim + dim
+end
+
+@inline function each_dof_idx(system::AbstractSystem)
+    return each_dof_idx(get_n_dim(system), 1:get_n_points(system))
+end
+
+@inline function each_loc_dof_idx(system::AbstractSystem)
+    return each_dof_idx(get_n_dim(system), 1:get_n_loc_points(system))
+end
+
+@inline function each_dof_idx(system::AbstractSystem, idxs::AbstractVector{<:Integer})
+    return each_dof_idx(get_n_dim(system), idxs)
+end
+
+#=
+This function generates a cartesian product of dof indices, dimensions, and point indices.
+In Julia it is very convenient, because elements in a 2-dimensional array can be
+accessed via:
+    A[dof]
+or via
+    A[dim, i]
+=#
+@inline function each_dof_idx(n_dim::Int, idxs::AbstractVector{<:Integer})
+    return ((get_dof(n_dim, dim, i), dim, i) for i in idxs, dim in 1:n_dim)
+end
+
+@inline function each_dof(system::AbstractSystem)
+    return each_dof(get_n_dim(system), 1:get_n_points(system))
 end
 
+@inline function each_loc_dof(system::AbstractSystem)
+    return each_dof(get_n_dim(system), 1:get_n_loc_points(system))
+end
+
+@inline function each_dof(system::AbstractSystem, idxs::AbstractVector{<:Integer})
+    return each_dof(get_n_dim(system), idxs)
+end
+
+@inline function each_dof(n_dim::Int, idxs::AbstractVector{<:Integer})
+    return (get_dof(n_dim, dim, i) for i in idxs, dim in 1:n_dim)
+end
+
+# Get the point index from a dof index.
+get_point(system::AbstractSystem, idx::Int) = get_point(get_n_dim(system), idx)
+get_point(n_dim::Int, idx::Int) = div(idx - 1, n_dim) + 1
+
+# Get the dimension index from a dof index.
+get_dim(system::AbstractSystem, idx::Int) = get_dim(get_n_dim(system), idx)
+get_dim(n_dim::Int, idx::Int) = mod(idx - 1, n_dim) + 1
+
 @inline function init_field_system(system, field)
     return nothing
 end

src/discretization/body.jl

@@ -91,9 +91,10 @@ struct Body{M<:AbstractMaterial,P<:AbstractPointParameters} <: AbstractBody{M}
     params_map::Vector{Int}
     single_dim_bcs::Vector{SingleDimBC}
     posdep_single_dim_bcs::Vector{PosDepSingleDimBC}
+    pos_single_dim_bcs::Vector{PosSingleDimBC}
+    data_bcs::Vector{DataBC}
     single_dim_ics::Vector{SingleDimIC}
     posdep_single_dim_ics::Vector{PosDepSingleDimIC}
-    data_bcs::Vector{DataBC}
     point_sets_precracks::Vector{PointSetsPreCrack}
 
     function Body(mat::M, position::AbstractMatrix, volume::AbstractVector) where {M}
@@ -109,14 +110,16 @@ struct Body{M<:AbstractMaterial,P<:AbstractPointParameters} <: AbstractBody{M}
 
         single_dim_bcs = Vector{SingleDimBC}()
         posdep_single_dim_bcs = Vector{PosDepSingleDimBC}()
+        pos_single_dim_bcs = Vector{PosSingleDimBC}()
         single_dim_ics = Vector{SingleDimIC}()
         posdep_single_dim_ics = Vector{PosDepSingleDimIC}()
-        v_bcs = Vector{DataBC}()
+        data_bcs = Vector{DataBC}()
         point_sets_precracks = Vector{PointSetsPreCrack}()
 
         new{M,P}(name, mat, n_points, position, volume, fail_permit, point_sets,
-                 point_params, params_map, single_dim_bcs, posdep_single_dim_bcs, v_bcs,
-                 single_dim_ics, posdep_single_dim_ics, point_sets_precracks)
+                 point_params, params_map, single_dim_bcs, posdep_single_dim_bcs,
+                 pos_single_dim_bcs, data_bcs, single_dim_ics, posdep_single_dim_ics,
+                 point_sets_precracks)
     end
 end
 
@@ -163,6 +166,10 @@ function Base.show(io::IO, ::MIME"text/plain", @nospecialize(body::AbstractBody)
             print(io, "\n    ")
             show(io, bc)
         end
+        for bc in body.pos_single_dim_bcs
+            print(io, "\n    ")
+            show(io, bc)
+        end
         for bc in body.data_bcs
             print(io, "\n    ")
             show(io, bc)
@@ -319,6 +326,11 @@ function log_msg_body(body::AbstractBody)
         settings = @sprintf("set `%s`, dimension %d", bc.point_set, bc.dim)
         msg *= msg_qty(descr, settings; indentation=4)
     end
+    for bc in body.pos_single_dim_bcs
+        descr = @sprintf("%s condition", field_to_name(bc.field))
+        settings = @sprintf("set `%s`, dimension %d", bc.point_set, bc.dim)
+        msg *= msg_qty(descr, settings; indentation=4)
+    end
     msg *= "  MATERIAL\n"
     msg *= log_material(body.mat; indentation=4)
     n_point_params = length(body.point_params)
@@ -372,8 +384,9 @@ end
 @inline function n_bcs(body::AbstractBody)
     n_sdbcs = length(body.single_dim_bcs)
     n_pdsdbcs = length(body.posdep_single_dim_bcs)
+    n_psdbcs = length(body.pos_single_dim_bcs)
     n_databcs = length(body.data_bcs)
-    return n_sdbcs + n_pdsdbcs + n_databcs
+    return n_sdbcs + n_pdsdbcs + n_psdbcs + n_databcs
 end
 @inline function n_ics(body::AbstractBody)
     n_sdics = length(body.single_dim_ics)

src/discretization/body_chunk.jl

@@ -34,10 +34,7 @@ struct BodyChunk{System<:AbstractSystem,
     mat::Material
     paramsetup::Params
     storage::Storage
-    psets::Dict{Symbol,Vector{Int}}
-    sdbcs::Vector{SingleDimBC}
-    pdsdbcs::Vector{PosDepSingleDimBC}
-    databcs::Vector{DataBC}
+    condhandler::ConditionHandler
     cells::Vector{MeshCell{VTKCellType,Tuple{Int64}}}
 end
 
@@ -48,13 +45,9 @@ function BodyChunk(body::AbstractBody, solver::AbstractTimeSolver, pd::PointDeco
     system = get_system(body, pd, chunk_id)
     paramsetup = get_paramsetup(body, system.chunk_handler, param_spec)
     storage = get_storage(mat, solver, system)
-    psets = localized_point_sets(body.point_sets, system.chunk_handler)
-    sdbcs = body.single_dim_bcs
-    pdsdbcs = body.posdep_single_dim_bcs
-    databcs = body.data_bcs
+    condhandler = ConditionHandler(body, system)
     cells = get_cells(get_n_loc_points(system))
-    chunk = BodyChunk(body_name, system, mat, paramsetup, storage, psets, sdbcs,
-                      pdsdbcs, databcs, cells)
+    chunk = BodyChunk(body_name, system, mat, paramsetup, storage, condhandler, cells)
     return chunk
 end
 
@@ -72,9 +65,20 @@ end
 
 @inline each_point_idx(chunk::AbstractBodyChunk) = each_point_idx(chunk.system)
 @inline each_point_idx_pair(chunk::AbstractBodyChunk) = each_point_idx_pair(chunk.system)
+@inline each_dof_idx(chunk::AbstractBodyChunk) = each_dof_idx(chunk.system)
+@inline each_loc_dof_idx(chunk::AbstractBodyChunk) = each_loc_dof_idx(chunk.system)
+@inline each_dof(chunk::AbstractBodyChunk) = each_dof(chunk.system)
+@inline each_loc_dof(chunk::AbstractBodyChunk) = each_loc_dof(chunk.system)
+@inline each_dim(chunk::AbstractBodyChunk) = each_dim(chunk.system)
 
-@inline n_loc_points(chunk::AbstractBodyChunk) = get_n_loc_points(chunk.system)
-@inline n_points(chunk::AbstractBodyChunk) = get_n_points(chunk.system)
+@inline get_n_loc_points(chunk::AbstractBodyChunk) = get_n_loc_points(chunk.system)
+@inline get_n_points(chunk::AbstractBodyChunk) = get_n_points(chunk.system)
+@inline get_n_dof(chunk::AbstractBodyChunk) = get_n_dof(chunk.system)
+@inline get_n_loc_dof(chunk::AbstractBodyChunk) = get_n_loc_dof(chunk.system)
+@inline get_n_dim(chunk::AbstractBodyChunk) = get_n_dim(chunk.system)
+
+@inline free_dofs(chunk::AbstractBodyChunk) = free_dofs(chunk.condhandler)
+@inline constrained_dofs(chunk::AbstractBodyChunk) = constrained_dofs(chunk.condhandler)
 
 function initialize!(::AbstractBodyChunk)
     return nothing

src/discretization/bond_associated_system.jl

@@ -14,12 +14,9 @@ end
 
 function BondAssociatedSystem(body::AbstractBody, pd::PointDecomposition, chunk_id::Int)
     check_bond_associated_system_compat(body.mat)
-    loc_points = pd.decomp[chunk_id]
-    bonds, n_neighbors = find_bonds(body, loc_points)
-    bond_ids = find_bond_ids(n_neighbors)
-    intersection_bond_ids = find_intersection_bond_ids(body, loc_points, bonds, bond_ids)
-    chunk_handler = get_chunk_handler(bonds, pd, chunk_id)
-    localize!(bonds, chunk_handler.localizer)
+    bonds, n_neighbors, bond_ids, chunk_handler = get_bond_data(body, pd, chunk_id)
+    intersection_bond_ids = find_intersection_bond_ids(body, chunk_handler.loc_points,
+                                                       bonds, bond_ids)
     position, volume = get_pos_and_vol_chunk(body, chunk_handler.point_ids)
     hood_volume = zeros(get_n_points(chunk_handler))
     ba_hood_volume = zeros(length(bonds))

src/discretization/bond_system.jl

@@ -54,10 +54,7 @@ end
 
 function BondSystem(body::AbstractBody, pd::PointDecomposition, chunk_id::Int)
     check_bond_system_compat(body.mat)
-    bonds, n_neighbors = find_bonds(body, pd.decomp[chunk_id])
-    bond_ids = find_bond_ids(n_neighbors)
-    chunk_handler = get_chunk_handler(bonds, pd, chunk_id)
-    localize!(bonds, chunk_handler.localizer)
+    bonds, n_neighbors, bond_ids, chunk_handler = get_bond_data(body, pd, chunk_id)
     position, volume = get_pos_and_vol_chunk(body, chunk_handler.point_ids)
     correction = get_correction(body.mat, chunk_handler.n_loc_points,
                                 length(chunk_handler.point_ids), length(bonds))
@@ -67,6 +64,16 @@ function BondSystem(body::AbstractBody, pd::PointDecomposition, chunk_id::Int)
     return system
 end
 
+function get_bond_data(body::AbstractBody, pd::PointDecomposition, chunk_id)
+    loc_points = pd.decomp[chunk_id]
+    bonds, n_neighbors = find_bonds(body, loc_points)
+    halo_points = find_halo_points(bonds, loc_points)
+    chunk_handler = ChunkHandler(pd, halo_points, chunk_id)
+    localize!(bonds, chunk_handler.localizer)
+    bond_ids = find_bond_ids(n_neighbors)
+    return bonds, n_neighbors, bond_ids, chunk_handler
+end
+
 function get_system(body::AbstractBody{Material}, pd::PointDecomposition,
                     chunk_id::Int) where {Material<:AbstractBondSystemMaterial}
     return BondSystem(body, pd, chunk_id)
@@ -86,7 +93,7 @@ function check_bond_system_compat(::AbstractBondSystemMaterial)
     return nothing
 end
 
-function find_bonds(body::AbstractBody, loc_points::UnitRange{Int})
+function find_bonds(body::AbstractBody, loc_points::AbstractVector{Int})
     δmax = maximum_horizon(body)
     nhs = GridNeighborhoodSearch{3}(search_radius=δmax, n_points=body.n_points)
     initialize_grid!(nhs, body.position)
@@ -125,15 +132,15 @@ end
 end
 
 function filter_bonds!(bonds::Vector{Bond}, n_neighbors::Vector{Int},
-                       loc_points::UnitRange{Int}, body::AbstractBody)
+                       loc_points::AbstractVector{Int}, body::AbstractBody)
     for crack in body.point_sets_precracks
         filter_bonds_by_crack!(bonds, n_neighbors, loc_points, crack, body)
     end
     return nothing
 end
 
 function filter_bonds_by_crack!(bonds::Vector{Bond}, n_neighbors::Vector{Int},
-                                loc_points::UnitRange{Int}, crack::PointSetsPreCrack,
+                                loc_points::AbstractVector{Int}, crack::PointSetsPreCrack,
                                 body::AbstractBody)
     filter_bonds(crack) || return nothing
     set_a, set_b = body.point_sets[crack.set_a], body.point_sets[crack.set_b]
@@ -174,17 +181,6 @@ function get_pos_and_vol_chunk(body::AbstractBody, point_ids::AbstractVector{<:I
     return position, volume
 end
 
-function get_chunk_handler(bonds::Vector{Bond}, pd::PointDecomposition, chunk_id::Int)
-    loc_points = pd.decomp[chunk_id]
-    n_loc_points = length(loc_points)
-    halo_points = find_halo_points(bonds, loc_points)
-    hidxs_by_src = sort_halo_by_src!(halo_points, pd.point_src, length(loc_points))
-    point_ids = vcat(loc_points, halo_points)
-    localizer = find_localizer(point_ids)
-    return ChunkHandler(n_loc_points, point_ids, loc_points, halo_points, hidxs_by_src,
-                        localizer)
-end
-
 function find_kernels(body::AbstractBody, chunk_handler::ChunkHandler, bonds::Vector{Bond},
                       bond_ids::Vector{UnitRange{Int}})
     hasproperty(body.mat, :kernel) || return Vector{Float64}()
@@ -208,7 +204,7 @@ end
     return system.kernels[bond_id]
 end
 
-function find_halo_points(bonds::Vector{Bond}, loc_points::UnitRange{Int})
+function find_halo_points(bonds::Vector{Bond}, loc_points::AbstractVector{Int})
     halo_points = Vector{Int}()
     for bond in bonds
         j = bond.neighbor
@@ -219,7 +215,7 @@ function find_halo_points(bonds::Vector{Bond}, loc_points::UnitRange{Int})
     return halo_points
 end
 
-@inline each_bond_idx(system::AbstractBondSystem, point_id::Int) = system.bond_ids[point_id]
+@inline each_bond_idx(system::AbstractBondSystem, i::Int) = system.bond_ids[i]
 
 function localize!(bonds::Vector{Bond}, localizer::Dict{Int,Int})
     for i in eachindex(bonds)

src/discretization/chunk_handler.jl

@@ -12,7 +12,8 @@ A type to handle a body chunk and its communication to other chunks.
 - `loc_points::UnitRange{Int}`: Indices of local points of the chunk.
 - `halo_points::Vector{Int}`: Indices of halo points of the chunk.
 - `hidxs_by_src::Dict{Int,UnitRange{Int}}`: Dict specifying the indices of halo Points
-    depending on the body chunk they belong to.
+    depending on the body chunk they belong to. So `body_chunk => indices`, with indices
+    being the indices of the halo points in `point_ids`.
 - `localizer::Dict{Int,Int}`: Localizes global indices to local indices in this chunk.
 """
 struct ChunkHandler <: AbstractChunkHandler
@@ -24,6 +25,17 @@ struct ChunkHandler <: AbstractChunkHandler
     localizer::Dict{Int,Int}
 end
 
+function ChunkHandler(pd::PointDecomposition, halo_points::Vector{Int}, chunk_id::Int)
+    loc_points = pd.decomp[chunk_id]
+    n_loc_points = length(loc_points)
+    hidxs_by_src = sort_halo_by_src!(halo_points, pd.point_src, length(loc_points))
+    point_ids = vcat(loc_points, halo_points)
+    localizer = find_localizer(point_ids)
+    chunk_handler = ChunkHandler(n_loc_points, point_ids, loc_points, halo_points,
+                                 hidxs_by_src, localizer)
+    return chunk_handler
+end
+
 for __field in fieldnames(ChunkHandler)
     local __funcname = Symbol("get_$(__field)")
     local __accessor_func = quote