From 30fa9d5d3971d81f436d5d6fb25c320020e20029 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 29 Mar 2023 15:07:29 +0200 Subject: [PATCH 01/16] First sweep, mostly for abstracting grid access via dof handlers. --- src/Dofs/ConstraintHandler.jl | 29 +++++++++++++++-------------- src/Dofs/DofHandler.jl | 12 ++++++------ src/Dofs/MixedDofHandler.jl | 2 +- src/Dofs/apply_analytical.jl | 14 +++++++------- src/Dofs/sparsity_pattern.jl | 4 ++-- src/Export/VTK.jl | 2 +- src/Grid/grid.jl | 15 ++++++++++++++- src/L2_projection.jl | 12 ++++++------ src/iterators.jl | 10 +++++----- 9 files changed, 57 insertions(+), 43 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 04ad21bb73..90afd95e55 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -242,7 +242,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet) if length(dbc.faces) == 0 @warn("adding Dirichlet Boundary Condition to set containing 0 entities") end - celltype = getcelltype(ch.dh.grid) + celltype = getcelltype(getgrid(ch.dh)) @assert isconcretetype(celltype) # Extract stuff for the field @@ -306,7 +306,7 @@ function add_prescribed_dof!(ch::ConstraintHandler, constrained_dof::Int, inhomo return ch end -function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) where {Index<:BoundaryIndex} +function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(getgrid(ch.dh)))) where {Index<:BoundaryIndex} local_face_dofs, local_face_dofs_offset = _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundarydof_indices(eltype(bcfaces))) copy!(dbc.local_face_dofs, local_face_dofs) @@ -352,13 +352,14 @@ function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, b return local_face_dofs, local_face_dofs_offset end -function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) - if interpolation !== default_interpolation(typeof(ch.dh.grid.cells[first(cellset)])) +function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(getgrid(ch.dh)))) + grid = getgrid(ch.dh) + if interpolation !== default_interpolation(typeof(getcells(grid, first(cellset)))) @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.") end ncomps = length(dbc.components) - nnodes = getnnodes(ch.dh.grid) + nnodes = getnnodes(grid) interpol_points = getnbasefunctions(interpolation) node_dofs = zeros(Int, ncomps, nnodes) visited = falses(nnodes) @@ -485,7 +486,7 @@ function _update!(inhomogeneities::Vector{Float64}, f::Function, nodes::Set{Int} dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T counter = 1 for (idx, nodenumber) in enumerate(nodeidxs) - x = dh.grid.nodes[nodenumber].x + x = getcoordinates(getnodes(getgrid(dh), nodenumber)) bc_value = f(x, time) @assert length(bc_value) == length(components) for v in bc_value @@ -513,13 +514,13 @@ function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler) for field in unique_fields nd = ndim(ch.dh, field) - data = zeros(Float64, nd, getnnodes(ch.dh.grid)) + data = zeros(Float64, nd, getnnodes(getgrid(ch.dh))) for dbc in ch.dbcs dbc.field_name != field && continue if eltype(dbc.faces) <: BoundaryIndex functype = boundaryfunction(eltype(dbc.faces)) for (cellidx, faceidx) in dbc.faces - for facenode in functype(ch.dh.grid.cells[cellidx])[faceidx] + for facenode in functype(getcells(getgrid(ch.dh), cellidx))[faceidx] for component in dbc.components data[component, facenode] = 1 end @@ -851,7 +852,7 @@ end function add!(ch::ConstraintHandler{<:MixedDofHandler}, dbc::Dirichlet) dbc_added = false for fh in ch.dh.fieldhandlers - if !isnothing(_find_field(fh, dbc.field_name)) && _in_cellset(ch.dh.grid, fh.cellset, dbc.faces; all=false) + if !isnothing(_find_field(fh, dbc.field_name)) && _in_cellset(getgrid(ch.dh), fh.cellset, dbc.faces; all=false) # Dofs in `dbc` not in `fh` will be removed, hence `dbc.faces` must be copied. # Recreating the `dbc` will create a copy of `dbc.faces`. # In this case, add! will warn, unless `warn_not_in_cellset=false` @@ -867,11 +868,11 @@ function add!(ch::ConstraintHandler{<:MixedDofHandler}, dbc::Dirichlet) end function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet; warn_not_in_cellset=true) - if warn_not_in_cellset && !(_in_cellset(ch.dh.grid, fh.cellset, dbc.faces; all=true)) + if warn_not_in_cellset && !(_in_cellset(getgrid(ch.dh), fh.cellset, dbc.faces; all=true)) @warn("You are trying to add a constraint a face/edge/node that is not in the cellset of the fieldhandler. This location will be skipped") end - celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset + celltype = getcelltype(getgrid(ch.dh), first(fh.cellset)) #Assume same celltype of all cells in fh.cellset # Extract stuff for the field field_idx = find_field(fh, dbc.field_name) @@ -986,7 +987,7 @@ function add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet) is_legacy = !isempty(pdbc.face_pairs) && isempty(pdbc.face_map) if is_legacy for (mset, iset) in pdbc.face_pairs - collect_periodic_faces!(pdbc.face_map, ch.dh.grid, mset, iset, identity) # TODO: Better transform + collect_periodic_faces!(pdbc.face_map, getgrid(ch.dh), mset, iset, identity) # TODO: Better transform end end field_idx = find_field(ch.dh, pdbc.field_name) @@ -1024,9 +1025,8 @@ end function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::Interpolation, field_dim::Int, offset::Int, is_legacy::Bool, rotation_matrix::Union{Matrix{T},Nothing}, ::Type{dof_map_t}, iterator_f::F) where {T, dof_map_t, F <: Function} - grid = ch.dh.grid + grid = getgrid(ch.dh) face_map = pdbc.face_map - Tx = typeof(first(ch.dh.grid.nodes).x) # Vec{D,T} # Indices of the local dofs for the faces local_face_dofs, local_face_dofs_offset = @@ -1098,6 +1098,7 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In "Dirichlet boundary condition on the relevant nodeset.", :PeriodicDirichlet) all_node_idxs = Set{Int}() + Tx = get_coordinate_type(getnodes(grid, getnodes(getcells(grid, first(face_map).mirror[1]), 1))) min_x = Tx(i -> typemax(eltype(Tx))) max_x = Tx(i -> typemin(eltype(Tx))) for facepair in face_map, faceidx in (facepair.mirror, facepair.image) diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 6c590d27f7..db7829e0ad 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -101,7 +101,7 @@ The field is added to all cells of the underlying grid. In case no interpolation given, the default interpolation of the grid's celltype is used. If the grid uses several celltypes, [`add!(dh::MixedDofHandler, fh::FieldHandler)`](@ref) must be used instead. """ -function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(dh.grid))) +function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(getgrid(dh)))) @assert !isclosed(dh) @assert !in(name, dh.field_names) push!(dh.field_names, name) @@ -111,7 +111,7 @@ function add!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_ end # Method for supporting dim=1 default -function add!(dh::DofHandler, name::Symbol, ip::Interpolation=default_interpolation(getcelltype(dh.grid))) +function add!(dh::DofHandler, name::Symbol, ip::Interpolation=default_interpolation(getcelltype(getgrid(dh)))) return add!(dh, name, 1, ip) end @@ -127,7 +127,7 @@ function __close!(dh::DofHandler{dim}) where {dim} # `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is # stored in vertexdict[v] # TODO: No need to allocate this vector for fields that don't have vertex dofs - vertexdicts = [zeros(Int, getnnodes(dh.grid)) for _ in 1:nfields(dh)] + vertexdicts = [zeros(Int, getnnodes(getgrid(dh))) for _ in 1:nfields(dh)] # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem # An edge is determined from two vertices, but we also need to store the direction @@ -161,7 +161,7 @@ function __close!(dh::DofHandler{dim}) where {dim} push!(dh.cell_dofs_offset, 1) # dofs for the first cell start at 1 # loop over all the cells, and distribute dofs for all the fields - for (ci, cell) in enumerate(getcells(dh.grid)) + for (ci, cell) in enumerate(getcells(getgrid(dh))) @debug println("cell #$ci") for field_idx in 1:nfields(dh) interpolation_info = interpolation_infos[field_idx] @@ -291,7 +291,7 @@ function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int) return global_dofs end -cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int) = cellcoords!(global_coords, dh.grid, i) +cellcoords!(global_coords::Vector{<:Vec}, dh::DofHandler, i::Int) = cellcoords!(global_coords, getgrid(dh), i) function reshape_to_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol) where T # make sure the field exists @@ -302,7 +302,7 @@ function reshape_to_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol) where field_dim = getfielddim(dh, field_idx) space_dim = field_dim == 2 ? 3 : field_dim - data = fill(zero(T), space_dim, getnnodes(dh.grid)) + data = fill(zero(T), space_dim, getnnodes(getgrid(dh))) reshape_field_data!(data, dh, u, offset, field_dim) diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index 6bc6cff4c8..fec55cbe3c 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -645,7 +645,7 @@ function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) return data end -function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=1:getncells(dh.grid)) where T +function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=1:getncells(getgrid(dh))) where T for cell in CellIterator(dh, cellset, UpdateFlags(; nodes=true, coords=false, dofs=true)) _celldofs = celldofs(cell) diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl index 5e488657eb..977983dfba 100644 --- a/src/Dofs/apply_analytical.jl +++ b/src/Dofs/apply_analytical.jl @@ -1,18 +1,18 @@ function _default_interpolations(dh::MixedDofHandler) fhs = dh.fieldhandlers - getcelltype(i) = typeof(getcells(dh.grid, first(fhs[i].cellset))) + getcelltype(i) = typeof(getcells(getgrid(dh), first(fhs[i].cellset))) ntuple(i -> default_interpolation(getcelltype(i)), length(fhs)) end function _default_interpolation(dh::DofHandler) - return default_interpolation(typeof(getcells(dh.grid, 1))) + return default_interpolation(typeof(getcells(getgrid(dh), 1))) end """ apply_analytical!( a::AbstractVector, dh::AbstractDofHandler, fieldname::Symbol, - f::Function, cellset=1:getncells(dh.grid)) + f::Function, cellset=1:getncells(getgrid(dh))) Apply a solution `f(x)` by modifying the values in the degree of freedom vector `a` pertaining to the field `fieldname` for all cells in `cellset`. @@ -32,7 +32,7 @@ This function can be used to apply initial conditions for time dependent problem """ function apply_analytical!( a::AbstractVector, dh::DofHandler, fieldname::Symbol, f::Function, - cellset = 1:getncells(dh.grid)) + cellset = 1:getncells(getgrid(dh))) fieldname ∉ getfieldnames(dh) && error("The fieldname $fieldname was not found in the dof handler") ip_geo = _default_interpolation(dh) @@ -45,7 +45,7 @@ end function apply_analytical!( a::AbstractVector, dh::MixedDofHandler, fieldname::Symbol, f::Function, - cellset = 1:getncells(dh.grid)) + cellset = 1:getncells(getgrid(dh))) fieldname ∉ getfieldnames(dh) && error("The fieldname $fieldname was not found in the dof handler") ip_geos = _default_interpolations(dh) @@ -65,7 +65,7 @@ function _apply_analytical!( a::Vector, dh::AbstractDofHandler, celldofinds, field_dim, ip_fun::Interpolation{dim,RefShape}, ip_geo::Interpolation, f::Function, cellset) where {dim, RefShape} - coords = getcoordinates(dh.grid, first(cellset)) + coords = getcoordinates(getgrid(dh), first(cellset)) ref_points = reference_coordinates(ip_fun) dummy_weights = zeros(length(ref_points)) qr = QuadratureRule{dim, RefShape}(dummy_weights, ref_points) @@ -77,7 +77,7 @@ function _apply_analytical!( length(f(first(coords))) == field_dim || error("length(f(x)) must be equal to dimension of the field ($field_dim)") for cellnr in cellset - getcoordinates!(coords, dh.grid, cellnr) + getcoordinates!(coords, getgrid(dh), cellnr) celldofs!(c_dofs, dh, cellnr) for (i, celldofind) in enumerate(celldofinds) f_dofs[i] = c_dofs[celldofind] diff --git a/src/Dofs/sparsity_pattern.jl b/src/Dofs/sparsity_pattern.jl index 21c8d5499f..7cbb1f25ff 100644 --- a/src/Dofs/sparsity_pattern.jl +++ b/src/Dofs/sparsity_pattern.jl @@ -113,7 +113,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint # of entries eliminated by constraints. max_buffer_length = ndofs(dh) # diagonal elements for (fhi, fh) in pairs(dh isa DofHandler ? (dh, ) : dh.fieldhandlers) - set = fh isa DofHandler ? (1:getncells(dh.grid)) : fh.cellset + set = fh isa DofHandler ? (1:getncells(getgrid(dh))) : fh.cellset n = ndofs_per_cell(dh, first(set)) # TODO: ndofs_per_cell(fh) entries_per_cell = if coupling === nothing sym ? div(n * (n + 1), 2) : n^2 @@ -130,7 +130,7 @@ function _create_sparsity_pattern(dh::AbstractDofHandler, ch#=::Union{Constraint for (fhi, fh) in pairs(dh isa DofHandler ? (dh, ) : dh.fieldhandlers) coupling === nothing || (coupling_fh = couplings[fhi]) - set = fh isa DofHandler ? (1:getncells(dh.grid)) : fh.cellset + set = fh isa DofHandler ? (1:getncells(getgrid(dh))) : fh.cellset n = ndofs_per_cell(dh, first(set)) # TODO: ndofs_per_cell(fh) resize!(global_dofs, n) @inbounds for element_id in set diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 9a1cf6de6e..2a178b6470 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -109,7 +109,7 @@ vtk_cellset(vtk::WriteVTK.DatasetFile, grid::AbstractGrid, cellset::String) = vtk_cellset(vtk, grid, [cellset]) function WriteVTK.vtk_grid(filename::AbstractString, dh::AbstractDofHandler; compress::Bool=true) - vtk_grid(filename, dh.grid; compress=compress) + vtk_grid(filename, getgrid(dh); compress=compress) end function WriteVTK.vtk_point_data(vtkfile, dh::AbstractDofHandler, u::Vector, suffix="") diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index ee862abab1..1f6d42546c 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -13,10 +13,23 @@ struct Node{dim,T} x::Vec{dim,T} end Node(x::NTuple{dim,T}) where {dim,T} = Node(Vec{dim,T}(x)) + +""" + getcoordinates(::Node) + +Get the value of the node coordinate. +""" getcoordinates(n::Node) = n.x """ - Ferrite.get_coordinate_eltype(::Node) + get_coordinate_type(::Node) + +Get the data type of the the node coordinate. +""" +get_coordinate_type(::Node{dim,T}) = Vec{dim,T} + +""" + get_coordinate_eltype(::Node) Get the data type of the components of the nodes coordinate. """ diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 353872e922..57ce847ade 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -104,7 +104,7 @@ end function Base.show(io::IO, ::MIME"text/plain", proj::L2Projector) println(io, typeof(proj)) - println(io, " projection on: ", length(proj.set), "/", getncells(proj.dh.grid), " cells in grid") + println(io, " projection on: ", length(proj.set), "/", getncells(getgrid(proj.dh)), " cells in grid") println(io, " function interpolation: ", proj.func_ip) println(io, " geometric interpolation: ", proj.geom_ip) end @@ -131,7 +131,7 @@ function _assemble_L2_matrix(fe_values, set, dh) celldofs!(cell_dofs, dh, cellnum) fill!(Me, 0) - Xe = getcoordinates(dh.grid, cellnum) + Xe = getcoordinates(getgrid(dh), cellnum) reinit!(fe_values, Xe) ## ∭( v ⋅ u )dΩ @@ -213,7 +213,7 @@ function project(proj::L2Projector, projected_vals = _project(vars, proj, fe_values, M, T)::Vector{T} if project_to_nodes # NOTE we may have more projected values than verticies in the mesh => not all values are returned - nnodes = getnnodes(proj.dh.grid) + nnodes = getnnodes(getgrid(proj.dh)) reordered_vals = fill(convert(T, NaN * zero(T)), nnodes) for node = 1:nnodes if (k = proj.node2dof_map[node]; k != 0) @@ -248,7 +248,7 @@ function _project(vars, proj::L2Projector, fe_values::Values, M::Integer, ::Type for (ic,cellnum) in enumerate(proj.set) celldofs!(cell_dofs, proj.dh, cellnum) fill!(fe, 0) - Xe = getcoordinates(proj.dh.grid, cellnum) + Xe = getcoordinates(getgrid(proj.dh), cellnum) cell_vars = vars[ic] reinit!(fe_values, Xe) @@ -279,7 +279,7 @@ end function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, proj::L2Projector, vals::Vector{T}, name::AbstractString) where T data = reshape_to_nodes(proj, vals) - @assert size(data, 2) == getnnodes(proj.dh.grid) + @assert size(data, 2) == getnnodes(getgrid(proj.dh)) vtk_point_data(vtk, data, name; component_names=component_names(T)) return vtk end @@ -295,7 +295,7 @@ function reshape_to_nodes(proj::L2Projector, vals::AbstractVector{S}) where {ord # can be any tensor field, however, the number of dofs should always match the length of vals @assert ndofs(dh) == length(vals) nout = S <: Vec{2} ? 3 : M # Pad 2D Vec to 3D - data = fill(T(NaN), nout, getnnodes(dh.grid)) + data = fill(T(NaN), nout, getnnodes(getgrid(dh))) for cell in CellIterator(dh, proj.set) _celldofs = celldofs(cell) @assert length(getnodes(cell)) == length(_celldofs) diff --git a/src/iterators.jl b/src/iterators.jl index a37a9f28c5..31c49d41bb 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -58,12 +58,12 @@ function CellCache(grid::Grid{dim,C,T}, flags::UpdateFlags=UpdateFlags()) where end function CellCache(dh::Union{DofHandler{dim},MixedDofHandler{dim}}, flags::UpdateFlags=UpdateFlags()) where {dim} - N = nnodes_per_cell(dh.grid) + N = nnodes_per_cell(getgrid(dh)) nodes = zeros(Int, N) - coords = zeros(Vec{dim, get_coordinate_eltype(dh.grid)}, N) + coords = zeros(Vec{dim, get_coordinate_eltype(getgrid(dh))}, N) n = ndofs_per_cell(dh) celldofs = zeros(Int, n) - return CellCache(flags, dh.grid, ScalarWrapper(-1), nodes, coords, dh, celldofs) + return CellCache(flags, getgrid(dh), ScalarWrapper(-1), nodes, coords, dh, celldofs) end # TODO: Can always resize and combine the two reinit! methods maybe? @@ -154,13 +154,13 @@ function CellIterator(gridordh::Union{Grid,AbstractDofHandler}, set::Union{IntegerCollection,Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) if set === nothing - grid = gridordh isa AbstractDofHandler ? gridordh.grid : gridordh + grid = gridordh isa AbstractDofHandler ? getgrid(gridordh) : gridordh set = 1:getncells(grid) end if gridordh isa MixedDofHandler # TODO: Since the CellCache is resizeable this is not really necessary to check # here, but might be useful to catch slow code paths? - _check_same_celltype(gridordh.grid, set) + _check_same_celltype(getgrid(gridordh), set) end return CellIterator(CellCache(gridordh, flags), set) end From ce09d344affb35af6f41b000b092a1870fee9b6a Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 29 Mar 2023 15:54:17 +0200 Subject: [PATCH 02/16] Fix oopsies. --- src/Dofs/ConstraintHandler.jl | 40 +++++++++++++++++------------------ src/Dofs/MixedDofHandler.jl | 7 ++++++ src/Grid/grid.jl | 2 +- 3 files changed, 28 insertions(+), 21 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 90afd95e55..70cb92d1f3 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -1106,14 +1106,14 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In nodes = faces(grid.cells[cellidx])[faceidx] union!(all_node_idxs, nodes) for n in nodes - x = grid.nodes[n].x + x = getcoordinates(getnodes(grid, n)) min_x = Tx(i -> min(min_x[i], x[i])) max_x = Tx(i -> max(max_x[i], x[i])) end end all_node_idxs_v = collect(all_node_idxs) points = construct_cornerish(min_x, max_x) - tree = KDTree(Tx[grid.nodes[i].x for i in all_node_idxs_v]) + tree = KDTree(Tx[getcoordinates(getnodes(grid, i)) for i in all_node_idxs_v]) idxs, _ = NearestNeighbors.nn(tree, points) corner_set = Set{Int}(all_node_idxs_v[i] for i in idxs) @@ -1387,12 +1387,12 @@ function __collect_periodic_faces_tree!(face_map::Vector{PeriodicFacePair}, grid if length(mset) != length(mset) error("different number of faces in mirror and image set") end - Tx = typeof(first(grid.nodes).x) + Tx = typeof(getcoordinates(first(getnodes(grid)))) mirror_mean_x = Tx[] for (c, f) in mset fn = faces(grid.cells[c])[f] - push!(mirror_mean_x, sum(grid.nodes[i].x for i in fn) / length(fn)) + push!(mirror_mean_x, sum(getcoordinates(getnodes(grid,i)) for i in fn) / length(fn)) end # Same dance for the image @@ -1400,7 +1400,7 @@ function __collect_periodic_faces_tree!(face_map::Vector{PeriodicFacePair}, grid for (c, f) in iset fn = faces(grid.cells[c])[f] # Apply transformation to all coordinates - push!(image_mean_x, sum(transformation(grid.nodes[i].x)::Tx for i in fn) / length(fn)) + push!(image_mean_x, sum(transformation(getcoordinates(getnodes(grid,i)))::Tx for i in fn) / length(fn)) end # Use KDTree to find closest face @@ -1477,16 +1477,16 @@ function __periodic_options(::T) where T <: Vec{3} end function __outward_normal(grid::Grid{2}, nodes, transformation::F=identity) where F <: Function - n1::Vec{2} = transformation(grid.nodes[nodes[1]].x) - n2::Vec{2} = transformation(grid.nodes[nodes[2]].x) + n1::Vec{2} = transformation(getcoordinates(getnodes(grid, nodes[1]))) + n2::Vec{2} = transformation(getcoordinates(getnodes(grid, nodes[2]))) n = Vec{2}((n2[2] - n1[2], - n2[1] + n1[1])) return n / norm(n) end function __outward_normal(grid::Grid{3}, nodes, transformation::F=identity) where F <: Function - n1::Vec{3} = transformation(grid.nodes[nodes[1]].x) - n2::Vec{3} = transformation(grid.nodes[nodes[2]].x) - n3::Vec{3} = transformation(grid.nodes[nodes[3]].x) + n1::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[1])))) + n2::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[2])))) + n3::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[3])))) n = (n3 - n2) × (n1 - n2) return n / norm(n) end @@ -1512,10 +1512,10 @@ function __check_periodic_faces(grid::Grid, fi::FaceIndex, fj::FaceIndex, known_ end # 2. Find the periodic direction using the vector between the midpoint of the faces - xmi = sum(grid.nodes[i].x for i in nodes_i) / length(nodes_i) - xmj = sum(grid.nodes[i].x for i in nodes_j) / length(nodes_j) + xmi = sum(getcoordinates(getnodes(grid, i)) for i in nodes_i) / length(nodes_i) + xmj = sum(getcoordinates(getnodes(grid, i)) for i in nodes_j) / length(nodes_j) xmij = xmj - xmi - h = 2 * norm(xmj - grid.nodes[nodes_j[1]].x) # Approximate element size + h = 2 * norm(xmj - getcoordinates(getnodes(grid, nodes_j[1]))) # Approximate element size TOLh = TOL * h found = false local len @@ -1531,11 +1531,11 @@ function __check_periodic_faces(grid::Grid, fi::FaceIndex, fj::FaceIndex, known_ # 3. Check that the first node of fj have a corresponding node in fi # In this method faces are mirrored (opposite normal vectors) so reverse the nodes nodes_i = circshift_tuple(reverse(nodes_i), 1) - xj = grid.nodes[nodes_j[1]].x + xj = getcoordinates(getnodes(grid, nodes_j[1])) node_rot = 0 found = false for i in eachindex(nodes_i) - xi = grid.nodes[nodes_i[i]].x + xi = getcoordinates(getnodes(grid, nodes_i[1])) xij = xj - xi if norm(xij - xmij) < TOLh found = true @@ -1547,8 +1547,8 @@ function __check_periodic_faces(grid::Grid, fi::FaceIndex, fj::FaceIndex, known_ # 4. Check the remaining nodes for the same criteria, now with known node_rot for j in 2:length(nodes_j) - xi = grid.nodes[nodes_i[mod1(j + node_rot, end)]].x - xj = grid.nodes[nodes_j[j]].x + xi = getcoordinates(getnodes(grid, nodes_i[mod1(j + node_rot, end)])) + xj = getcoordinates(getnodes(grid, nodes_j[j])) xij = xj - xi if norm(xij - xmij) >= TOLh return nothing @@ -1594,14 +1594,14 @@ function __check_periodic_faces_f(grid::Grid, fi::FaceIndex, fj::FaceIndex, xmi, # 2. Compute the relative rotation xmij = xmj - xmi - h = 2 * norm(xmj - grid.nodes[nodes_j[1]].x) # Approximate element size + h = 2 * norm(xmj - getcoordinates(getnodes(grid, nodes_j[1]))) # Approximate element size TOLh = TOL * h nodes_i = mirror ? circshift_tuple(reverse(nodes_i), 1) : nodes_i # reverse if necessary - xj = transformation(grid.nodes[nodes_j[1]].x) + xj = transformation(getcoordinates(getnodes(grid, nodes_j[1]))) node_rot = 0 found = false for i in eachindex(nodes_i) - xi = grid.nodes[nodes_i[i]].x + xi = getcoordinates(getnodes(grid, nodes_i[i])) xij = xj - xi if norm(xij - xmij) < TOLh found = true diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index fec55cbe3c..57ac233a67 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -74,6 +74,13 @@ end isclosed(dh::AbstractDofHandler) = dh.closed[] +""" + getgrid(dh::AbstractDofHandler) + +Return the grid on which the dof handler acts. +""" +getgrid(dh::AbstractDofHandler) = dh.grid + """ ndofs(dh::AbstractDofHandler) diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 1f6d42546c..0350cb4907 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -26,7 +26,7 @@ getcoordinates(n::Node) = n.x Get the data type of the the node coordinate. """ -get_coordinate_type(::Node{dim,T}) = Vec{dim,T} +get_coordinate_type(::Node{dim,T}) where {dim,T} = Vec{dim,T} """ get_coordinate_eltype(::Node) From 28350121a744ec9d55299434bed5a3cf7f720fe7 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 29 Mar 2023 16:00:51 +0200 Subject: [PATCH 03/16] More cleanup. --- src/Dofs/ConstraintHandler.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 70cb92d1f3..3bd98422ce 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -1387,7 +1387,7 @@ function __collect_periodic_faces_tree!(face_map::Vector{PeriodicFacePair}, grid if length(mset) != length(mset) error("different number of faces in mirror and image set") end - Tx = typeof(getcoordinates(first(getnodes(grid)))) + Tx = get_coordinate_type(first(getnodes(grid))) mirror_mean_x = Tx[] for (c, f) in mset @@ -1484,9 +1484,9 @@ function __outward_normal(grid::Grid{2}, nodes, transformation::F=identity) wher end function __outward_normal(grid::Grid{3}, nodes, transformation::F=identity) where F <: Function - n1::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[1])))) - n2::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[2])))) - n3::Vec{3} = transformation(transformation(getcoordinates(getnodes(grid, nodes[3])))) + n1::Vec{3} = transformation(getcoordinates(getnodes(grid, nodes[1]))) + n2::Vec{3} = transformation(getcoordinates(getnodes(grid, nodes[2]))) + n3::Vec{3} = transformation(getcoordinates(getnodes(grid, nodes[3]))) n = (n3 - n2) × (n1 - n2) return n / norm(n) end @@ -1513,7 +1513,7 @@ function __check_periodic_faces(grid::Grid, fi::FaceIndex, fj::FaceIndex, known_ # 2. Find the periodic direction using the vector between the midpoint of the faces xmi = sum(getcoordinates(getnodes(grid, i)) for i in nodes_i) / length(nodes_i) - xmj = sum(getcoordinates(getnodes(grid, i)) for i in nodes_j) / length(nodes_j) + xmj = sum(getcoordinates(getnodes(grid, j)) for j in nodes_j) / length(nodes_j) xmij = xmj - xmi h = 2 * norm(xmj - getcoordinates(getnodes(grid, nodes_j[1]))) # Approximate element size TOLh = TOL * h @@ -1535,7 +1535,7 @@ function __check_periodic_faces(grid::Grid, fi::FaceIndex, fj::FaceIndex, known_ node_rot = 0 found = false for i in eachindex(nodes_i) - xi = getcoordinates(getnodes(grid, nodes_i[1])) + xi = getcoordinates(getnodes(grid, nodes_i[i])) xij = xj - xi if norm(xij - xmij) < TOLh found = true From fc3b69e10808e20b7a5b7c1f99ee3d6714e3dea7 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 29 Mar 2023 16:50:30 +0200 Subject: [PATCH 04/16] I think keeping the assumption that grids come with nodal spaces is reasonable - rolling back one change. --- src/Dofs/ConstraintHandler.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 3bd98422ce..944b2074c6 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -1098,7 +1098,7 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In "Dirichlet boundary condition on the relevant nodeset.", :PeriodicDirichlet) all_node_idxs = Set{Int}() - Tx = get_coordinate_type(getnodes(grid, getnodes(getcells(grid, first(face_map).mirror[1]), 1))) + Tx = get_coordinate_type(getnodes(grid, getcells(grid, 1).nodes[1])) min_x = Tx(i -> typemax(eltype(Tx))) max_x = Tx(i -> typemin(eltype(Tx))) for facepair in face_map, faceidx in (facepair.mirror, facepair.image) From f18d86b8421c1357add792c94a543c563d7924c7 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 29 Mar 2023 22:35:22 +0200 Subject: [PATCH 05/16] Get coordinate type directly from grid. --- src/Dofs/ConstraintHandler.jl | 4 ++-- src/Grid/grid.jl | 7 +++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 944b2074c6..771cc26c06 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -1098,7 +1098,7 @@ function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::In "Dirichlet boundary condition on the relevant nodeset.", :PeriodicDirichlet) all_node_idxs = Set{Int}() - Tx = get_coordinate_type(getnodes(grid, getcells(grid, 1).nodes[1])) + Tx = get_coordinate_type(grid) min_x = Tx(i -> typemax(eltype(Tx))) max_x = Tx(i -> typemin(eltype(Tx))) for facepair in face_map, faceidx in (facepair.mirror, facepair.image) @@ -1387,7 +1387,7 @@ function __collect_periodic_faces_tree!(face_map::Vector{PeriodicFacePair}, grid if length(mset) != length(mset) error("different number of faces in mirror and image set") end - Tx = get_coordinate_type(first(getnodes(grid))) + Tx = get_coordinate_type(grid) mirror_mean_x = Tx[] for (c, f) in mset diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index 0350cb4907..d04a705fc9 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -410,6 +410,13 @@ end ########################## # Grid utility functions # ########################## +""" + get_coordinate_type(::AbstractGrid) + +Get the datatype for a single point in the grid. +""" +get_coordinate_type(grid::Grid{dim,C,T}) where {dim,C,T} = Vec{dim,T} # Node is baked into the mesh type. + """ getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, cellidx::CellIndex, include_self=false) getneighborhood(top::ExclusiveTopology, grid::AbstractGrid, faceidx::FaceIndex, include_self=false) From aa4cc42573d1f8f552b8ff0f3cfeadf49718b6d0 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Thu, 30 Mar 2023 11:16:12 +0200 Subject: [PATCH 06/16] Apply suggestions from code review Co-authored-by: Knut Andreas Meyer --- src/Dofs/ConstraintHandler.jl | 2 +- src/Dofs/apply_analytical.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 771cc26c06..6ce543669c 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -354,7 +354,7 @@ end function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(getgrid(ch.dh)))) grid = getgrid(ch.dh) - if interpolation !== default_interpolation(typeof(getcells(grid, first(cellset)))) + if interpolation !== default_interpolation(getcelltype(grid, first(cellset))) @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.") end diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl index 977983dfba..cdde179f28 100644 --- a/src/Dofs/apply_analytical.jl +++ b/src/Dofs/apply_analytical.jl @@ -1,12 +1,12 @@ function _default_interpolations(dh::MixedDofHandler) fhs = dh.fieldhandlers - getcelltype(i) = typeof(getcells(getgrid(dh), first(fhs[i].cellset))) + getcelltype(i) = getcelltype(getgrid(dh), first(fhs[i].cellset)) ntuple(i -> default_interpolation(getcelltype(i)), length(fhs)) end function _default_interpolation(dh::DofHandler) - return default_interpolation(typeof(getcells(getgrid(dh), 1))) + return default_interpolation(getcelltype(getgrid(dh), 1)) end """ From 4cdc6979a35630df94171d50b9f2f7f2935be398 Mon Sep 17 00:00:00 2001 From: termi-official Date: Thu, 30 Mar 2023 12:43:54 +0200 Subject: [PATCH 07/16] Fix oopsie and propagate interface in remaining code (hopefully). --- src/Dofs/MixedDofHandler.jl | 29 ++++++++++++++++++----------- src/Dofs/apply_analytical.jl | 5 ++--- src/Export/VTK.jl | 2 +- src/Grid/grid.jl | 2 +- src/L2_projection.jl | 4 ++-- src/iterators.jl | 6 +++--- 6 files changed, 27 insertions(+), 21 deletions(-) diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl index 57ac233a67..903d726453 100644 --- a/src/Dofs/MixedDofHandler.jl +++ b/src/Dofs/MixedDofHandler.jl @@ -96,7 +96,14 @@ Return the number of degrees of freedom for the cell with index `cell`. See also [`ndofs`](@ref). """ ndofs_per_cell(dh::MixedDofHandler, cell::Int=1) = dh.ndofs_per_cell[cell] -nnodes_per_cell(dh::MixedDofHandler, cell::Int=1) = nnodes_per_cell(dh.grid, cell) # TODO: deprecate, shouldn't belong to MixedDofHandler any longer +nnodes_per_cell(dh::MixedDofHandler, cell::Int=1) = nnodes_per_cell(getgrid(dh), cell) # TODO: deprecate, shouldn't belong to MixedDofHandler any longer + +""" + getcelltype(dh::MixedDofHandler, fh::FieldHandler) + +Shortcut to get the cell type of the field handler. +""" +getcelltype(dh::MixedDofHandler, fh::FieldHandler) = getcelltype(getgrid(dh), first(fh.cellset)) """ celldofs!(global_dofs::Vector{Int}, dh::AbstractDofHandler, i::Int) @@ -125,11 +132,11 @@ end #TODO: perspectively remove in favor of `getcoordinates!(global_coords, grid, i)`? function cellcoords!(global_coords::Vector{Vec{dim,T}}, dh::MixedDofHandler, i::Union{Int, <:AbstractCell}) where {dim,T} - cellcoords!(global_coords, dh.grid, i) + cellcoords!(global_coords, getgrid(dh), i) end function cellnodes!(global_nodes::Vector{Int}, dh::MixedDofHandler, i::Union{Int, <:AbstractCell}) - cellnodes!(global_nodes, dh.grid, i) + cellnodes!(global_nodes, getgrid(dh), i) end """ @@ -174,11 +181,11 @@ Add all fields of the [`FieldHandler`](@ref) `fh` to `dh`. function add!(dh::MixedDofHandler, fh::FieldHandler) # TODO: perhaps check that a field with the same name is the same field? @assert !isclosed(dh) - _check_same_celltype(dh.grid, collect(fh.cellset)) + _check_same_celltype(getgrid(dh), collect(fh.cellset)) _check_cellset_intersections(dh, fh) # the field interpolations should have the same refshape as the cells they are applied to # extract the celltype from the first cell as the celltypes are all equal - cell_type = typeof(dh.grid.cells[first(fh.cellset)]) + cell_type = getcelltype(getgrid(dh), first(fh.cellset)) refshape_cellset = getrefshape(default_interpolation(cell_type)) for field_idx in eachindex(fh.fields) refshape = getrefshape(getfieldinterpolation(fh, field_idx)) @@ -196,7 +203,7 @@ function _check_cellset_intersections(dh::MixedDofHandler, fh::FieldHandler) end function add!(dh::MixedDofHandler, name::Symbol, dim::Int) - celltype = getcelltype(dh.grid) + celltype = getcelltype(getgrid(dh)) isconcretetype(celltype) || error("If you have more than one celltype in Grid, you must use add!(dh::MixedDofHandler, fh::FieldHandler)") add!(dh, name, dim, default_interpolation(celltype)) end @@ -204,11 +211,11 @@ end function add!(dh::MixedDofHandler, name::Symbol, dim::Int, ip::Interpolation) @assert !isclosed(dh) - celltype = getcelltype(dh.grid) + celltype = getcelltype(getgrid(dh)) @assert isconcretetype(celltype) if length(dh.fieldhandlers) == 0 - cellset = Set(1:getncells(dh.grid)) + cellset = Set(1:getncells(getgrid(dh))) push!(dh.fieldhandlers, FieldHandler(Field[], cellset)) elseif length(dh.fieldhandlers) > 1 error("If you have more than one FieldHandler, you must specify field") @@ -252,7 +259,7 @@ function __close!(dh::MixedDofHandler{dim}) where {dim} # `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is # stored in vertexdict[v] # TODO: No need to allocate this vector for fields that don't have vertex dofs - vertexdicts = [zeros(Int, getnnodes(dh.grid)) for _ in 1:numfields] + vertexdicts = [zeros(Int, getnnodes(getgrid(dh))) for _ in 1:numfields] # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem. # An edge is uniquely determined by two vertices, but we also need to store the @@ -308,7 +315,7 @@ function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fiel for ci in cellnumbers @debug "Creating dofs for cell #$ci" - cell = dh.grid.cells[ci] + cell = getcells(getgrid(dh), ci) len_cell_dofs = length(dh.cell_dofs) dh.cell_dofs_offset[ci] = len_cell_dofs + 1 @@ -639,7 +646,7 @@ function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) field_dim = getfielddim(dh, fieldname) space_dim = field_dim == 2 ? 3 : field_dim - data = fill(T(NaN), space_dim, getnnodes(dh.grid)) # set default value + data = fill(T(NaN), space_dim, getnnodes(getgrid(dh))) # set default value for fh in dh.fieldhandlers # check if this fh contains this field, otherwise continue to the next diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl index cdde179f28..51e28685bd 100644 --- a/src/Dofs/apply_analytical.jl +++ b/src/Dofs/apply_analytical.jl @@ -1,12 +1,11 @@ function _default_interpolations(dh::MixedDofHandler) fhs = dh.fieldhandlers - getcelltype(i) = getcelltype(getgrid(dh), first(fhs[i].cellset)) - ntuple(i -> default_interpolation(getcelltype(i)), length(fhs)) + ntuple(i -> default_interpolation(getcelltype(dh, fhs[i])), length(fhs)) end function _default_interpolation(dh::DofHandler) - return default_interpolation(getcelltype(getgrid(dh), 1)) + return default_interpolation(getcelltype(getgrid(dh))) end """ diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index 2a178b6470..d0a6156679 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -26,7 +26,7 @@ which data can be appended to, see `vtk_point_data` and `vtk_cell_data`. """ function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; compress::Bool=true) where {dim,C,T} cls = MeshCell[] - for cell in grid.cells + for cell in getcells(grid) celltype = Ferrite.cell_to_vtkcell(typeof(cell)) push!(cls, MeshCell(celltype, nodes_to_vtkorder(cell))) end diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index d04a705fc9..0f921132e2 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -600,7 +600,7 @@ Returns all vertex sets of the grid. """ @inline getvertexsets(grid::AbstractGrid) = grid.vertexsets -n_faces_per_cell(grid::Grid) = nfaces(eltype(grid.cells)) +n_faces_per_cell(grid::Grid) = nfaces(getcelltype(grid)) """ function compute_vertex_values(grid::AbstractGrid, f::Function) diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 57ce847ade..539969ac13 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -35,7 +35,7 @@ end function L2Projector(qr::QuadratureRule, func_ip::Interpolation, grid::Ferrite.AbstractGrid, set=1:getncells(grid), qr_mass::QuadratureRule=_mass_qr(func_ip), - geom_ip::Interpolation = default_interpolation(typeof(grid.cells[first(set)]))) + geom_ip::Interpolation = default_interpolation(getcelltype(grid, first(set)))) Base.depwarn("L2Projector(qr, func_ip, grid) is deprecated, " * "use L2Projector(func_ip, grid) instead.", :L2Projector) return L2Projector(func_ip, grid; qr_lhs=qr_mass, set=set, geom_ip=geom_ip, qr_rhs=qr) @@ -69,7 +69,7 @@ function L2Projector( grid::AbstractGrid; qr_lhs::QuadratureRule = _mass_qr(func_ip), set = 1:getncells(grid), - geom_ip::Interpolation = default_interpolation(typeof(grid.cells[first(set)])), + geom_ip::Interpolation = default_interpolation(getcelltype(grid, first(set))), qr_rhs::Union{QuadratureRule,Nothing}=nothing, # deprecated ) diff --git a/src/iterators.jl b/src/iterators.jl index 31c49d41bb..cdf4c24192 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -113,7 +113,7 @@ cellid(cc::CellCache) = cc.cellid[] celldofs!(v::Vector, cc::CellCache) = copyto!(v, cc.dofs) # celldofs!(v, cc.dh, cc.cellid[]) # TODO: These should really be replaced with something better... -nfaces(cc::CellCache) = nfaces(eltype(cc.grid.cells)) +nfaces(cc::CellCache) = nfaces(getcelltype(cc.grid)) onboundary(cc::CellCache, face::Int) = cc.grid.boundary_matrix[face, cc.cellid[]] ################## @@ -183,8 +183,8 @@ Base.length(ci::CellIterator) = length(ci.set) function _check_same_celltype(grid::AbstractGrid, cellset) - celltype = typeof(grid.cells[first(cellset)]) - if !all(typeof(grid.cells[i]) == celltype for i in cellset) + celltype = getcelltype(grid, first(cellset)) + if !all(getcelltype(grid, i) == celltype for i in cellset) error("The cells in the cellset are not all of the same celltype.") end end From ad8941a4d33183ef8ab74cbd1cf9318fdf700d6d Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 7 Jun 2023 17:33:46 +0200 Subject: [PATCH 08/16] Export getgrid. --- src/exports.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/exports.jl b/src/exports.jl index 627b3f32c6..aaf47e869d 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -119,6 +119,7 @@ export Field, evaluate_at_grid_nodes, apply_analytical!, + getgrid, # Constraints ConstraintHandler, From cc509df7328b907ec462e013b82c3323e63282f0 Mon Sep 17 00:00:00 2001 From: termi-official Date: Wed, 7 Jun 2023 17:40:52 +0200 Subject: [PATCH 09/16] More exports. --- src/exports.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/exports.jl b/src/exports.jl index aaf47e869d..5293d122b5 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -120,6 +120,7 @@ export evaluate_at_grid_nodes, apply_analytical!, getgrid, + getdim, # Constraints ConstraintHandler, From 5914eab9125d37dc2f94853e45eafe5eb707fb2b Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Fri, 9 Jun 2023 14:18:42 +0200 Subject: [PATCH 10/16] Derp. --- src/Dofs/MixedDofHandler.jl | 680 ------------------------------------ 1 file changed, 680 deletions(-) delete mode 100644 src/Dofs/MixedDofHandler.jl diff --git a/src/Dofs/MixedDofHandler.jl b/src/Dofs/MixedDofHandler.jl deleted file mode 100644 index 903d726453..0000000000 --- a/src/Dofs/MixedDofHandler.jl +++ /dev/null @@ -1,680 +0,0 @@ -abstract type AbstractDofHandler end - -""" - Field(name::Symbol, interpolation::Interpolation, dim::Int) - -Construct `dim`-dimensional `Field` called `name` which is approximated by `interpolation`. - -The interpolation is used for distributing the degrees of freedom. -""" -struct Field - name::Symbol - interpolation::Interpolation - dim::Int -end - -""" - FieldHandler(fields::Vector{Field}, cellset::Set{Int}) - -Construct a `FieldHandler` based on an array of [`Field`](@ref)s and assigns it a set of cells. - -A `FieldHandler` must fullfill the following requirements: -- All [`Cell`](@ref)s in `cellset` are of the same type. -- Each field only uses a single interpolation on the `cellset`. -- Each cell belongs only to a single `FieldHandler`, i.e. all fields on a cell must be added within the same `FieldHandler`. - -Notice that a `FieldHandler` can hold several fields. -""" -mutable struct FieldHandler - fields::Vector{Field} - cellset::Set{Int} -end - -""" - MixedDofHandler(grid::Grid) - -Construct a `MixedDofHandler` based on `grid`. Supports: -- `Grid`s with or without concrete element type (E.g. "mixed" grids with several different element types.) -- One or several fields, which can live on the whole domain or on subsets of the `Grid`. -""" -struct MixedDofHandler{dim,T,G<:AbstractGrid{dim}} <: AbstractDofHandler - fieldhandlers::Vector{FieldHandler} - field_names::Vector{Symbol} - # Dofs for cell i are stored in cell_dofs[cell_dofs_offset[i]:(cell_dofs_offset[i]+length[i]-1)]. - # Note that explicitly keeping track of ndofs_per_cell is necessary since dofs are *not* - # distributed in cell order like for the DofHandler (where the length can be determined - # by cell_dofs_offset[i+1]-cell_dofs_offset[i]). - # TODO: ndofs_per_cell should probably be replaced by ndofs_per_fieldhandler, since all - # cells in a FieldHandler have the same number of dofs. - cell_dofs::Vector{Int} - cell_dofs_offset::Vector{Int} - ndofs_per_cell::Vector{Int} - closed::ScalarWrapper{Bool} - grid::G - ndofs::ScalarWrapper{Int} -end - -function MixedDofHandler(grid::Grid{dim,C,T}) where {dim,C,T} - ncells = getncells(grid) - MixedDofHandler{dim,T,typeof(grid)}(FieldHandler[], Symbol[], Int[], zeros(Int, ncells), zeros(Int, ncells), ScalarWrapper(false), grid, ScalarWrapper(-1)) -end - -function Base.show(io::IO, ::MIME"text/plain", dh::MixedDofHandler) - println(io, typeof(dh)) - println(io, " Fields:") - for fieldname in getfieldnames(dh) - println(io, " ", repr(fieldname), ", dim: ", getfielddim(dh, fieldname)) - end - if !isclosed(dh) - print(io, " Not closed!") - else - print(io, " Total dofs: ", ndofs(dh)) - end -end - -isclosed(dh::AbstractDofHandler) = dh.closed[] - -""" - getgrid(dh::AbstractDofHandler) - -Return the grid on which the dof handler acts. -""" -getgrid(dh::AbstractDofHandler) = dh.grid - -""" - ndofs(dh::AbstractDofHandler) - -Return the number of degrees of freedom in `dh` -""" -ndofs(dh::AbstractDofHandler) = dh.ndofs[] - -""" - ndofs_per_cell(dh::AbstractDofHandler[, cell::Int=1]) - -Return the number of degrees of freedom for the cell with index `cell`. - -See also [`ndofs`](@ref). -""" -ndofs_per_cell(dh::MixedDofHandler, cell::Int=1) = dh.ndofs_per_cell[cell] -nnodes_per_cell(dh::MixedDofHandler, cell::Int=1) = nnodes_per_cell(getgrid(dh), cell) # TODO: deprecate, shouldn't belong to MixedDofHandler any longer - -""" - getcelltype(dh::MixedDofHandler, fh::FieldHandler) - -Shortcut to get the cell type of the field handler. -""" -getcelltype(dh::MixedDofHandler, fh::FieldHandler) = getcelltype(getgrid(dh), first(fh.cellset)) - -""" - celldofs!(global_dofs::Vector{Int}, dh::AbstractDofHandler, i::Int) - -Store the degrees of freedom that belong to cell `i` in `global_dofs`. - -See also [`celldofs`](@ref). -""" -function celldofs!(global_dofs::Vector{Int}, dh::MixedDofHandler, i::Int) - @assert isclosed(dh) - @assert length(global_dofs) == ndofs_per_cell(dh, i) - unsafe_copyto!(global_dofs, 1, dh.cell_dofs, dh.cell_dofs_offset[i], length(global_dofs)) - return global_dofs -end - -""" - celldofs(dh::AbstractDofHandler, i::Int) - -Return a vector with the degrees of freedom that belong to cell `i`. - -See also [`celldofs!`](@ref). -""" -function celldofs(dh::AbstractDofHandler, i::Int) - return celldofs!(zeros(Int, ndofs_per_cell(dh, i)), dh, i) -end - -#TODO: perspectively remove in favor of `getcoordinates!(global_coords, grid, i)`? -function cellcoords!(global_coords::Vector{Vec{dim,T}}, dh::MixedDofHandler, i::Union{Int, <:AbstractCell}) where {dim,T} - cellcoords!(global_coords, getgrid(dh), i) -end - -function cellnodes!(global_nodes::Vector{Int}, dh::MixedDofHandler, i::Union{Int, <:AbstractCell}) - cellnodes!(global_nodes, getgrid(dh), i) -end - -""" - getfieldnames(dh::MixedDofHandler) - -Return a vector with the names of all fields. Can be used as an iterable over all the fields -in the problem. -""" -getfieldnames(dh::MixedDofHandler) = dh.field_names - -getfielddim(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].dim -getfielddim(fh::FieldHandler, field_name::Symbol) = getfielddim(fh, find_field(fh, field_name)) - -""" - getfielddim(dh::MixedDofHandler, field_idxs::NTuple{2,Int}) - getfielddim(dh::MixedDofHandler, field_name::Symbol) - getfielddim(dh::FieldHandler, field_idx::Int) - getfielddim(dh::FieldHandler, field_name::Symbol) - -Return the dimension of a given field. The field can be specified by its index (see -[`find_field`](@ref)) or its name. -""" -function getfielddim(dh::MixedDofHandler, field_idxs::NTuple{2, Int}) - fh_idx, field_idx = field_idxs - fielddim = getfielddim(dh.fieldhandlers[fh_idx], field_idx) - return fielddim -end -getfielddim(dh::MixedDofHandler, name::Symbol) = getfielddim(dh, find_field(dh, name)) - -""" - nfields(dh::MixedDofHandler) - -Returns the number of unique fields defined. -""" -nfields(dh::MixedDofHandler) = length(getfieldnames(dh)) - -""" - add!(dh::MixedDofHandler, fh::FieldHandler) - -Add all fields of the [`FieldHandler`](@ref) `fh` to `dh`. -""" -function add!(dh::MixedDofHandler, fh::FieldHandler) - # TODO: perhaps check that a field with the same name is the same field? - @assert !isclosed(dh) - _check_same_celltype(getgrid(dh), collect(fh.cellset)) - _check_cellset_intersections(dh, fh) - # the field interpolations should have the same refshape as the cells they are applied to - # extract the celltype from the first cell as the celltypes are all equal - cell_type = getcelltype(getgrid(dh), first(fh.cellset)) - refshape_cellset = getrefshape(default_interpolation(cell_type)) - for field_idx in eachindex(fh.fields) - refshape = getrefshape(getfieldinterpolation(fh, field_idx)) - refshape_cellset == refshape || error("The RefShapes of the fieldhandlers interpolations must correspond to the RefShape of the cells it is applied to.") - end - - push!(dh.fieldhandlers, fh) - return dh -end - -function _check_cellset_intersections(dh::MixedDofHandler, fh::FieldHandler) - for _fh in dh.fieldhandlers - isdisjoint(_fh.cellset, fh.cellset) || error("Each cell can only belong to a single FieldHandler.") - end -end - -function add!(dh::MixedDofHandler, name::Symbol, dim::Int) - celltype = getcelltype(getgrid(dh)) - isconcretetype(celltype) || error("If you have more than one celltype in Grid, you must use add!(dh::MixedDofHandler, fh::FieldHandler)") - add!(dh, name, dim, default_interpolation(celltype)) -end - -function add!(dh::MixedDofHandler, name::Symbol, dim::Int, ip::Interpolation) - @assert !isclosed(dh) - - celltype = getcelltype(getgrid(dh)) - @assert isconcretetype(celltype) - - if length(dh.fieldhandlers) == 0 - cellset = Set(1:getncells(getgrid(dh))) - push!(dh.fieldhandlers, FieldHandler(Field[], cellset)) - elseif length(dh.fieldhandlers) > 1 - error("If you have more than one FieldHandler, you must specify field") - end - fh = first(dh.fieldhandlers) - - field = Field(name,ip,dim) - - push!(fh.fields, field) - - return dh -end - -""" - close!(dh::AbstractDofHandler) - -Closes `dh` and creates degrees of freedom for each cell. - -If there are several fields, the dofs are added in the following order: -For a `MixedDofHandler`, go through each `FieldHandler` in the order they were added. -For each field in the `FieldHandler` or in the `DofHandler` (again, in the order the fields were added), -create dofs for the cell. -This means that dofs on a particular cell, the dofs will be numbered according to the fields; -first dofs for field 1, then field 2, etc. -""" -function close!(dh::MixedDofHandler) - dh, _, _, _ = __close!(dh) - return dh -end - -function __close!(dh::MixedDofHandler{dim}) where {dim} - @assert !isclosed(dh) - - # Collect the global field names - empty!(dh.field_names) - for fh in dh.fieldhandlers, f in fh.fields - f.name in dh.field_names || push!(dh.field_names, f.name) - end - numfields = length(dh.field_names) - - # `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is - # stored in vertexdict[v] - # TODO: No need to allocate this vector for fields that don't have vertex dofs - vertexdicts = [zeros(Int, getnnodes(getgrid(dh))) for _ in 1:numfields] - - # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem. - # An edge is uniquely determined by two vertices, but we also need to store the - # direction of the first edge we encounter and add dofs too. When we encounter the same - # edge the next time we check if the direction is the same, otherwise we reuse the dofs - # in the reverse order. - edgedicts = [Dict{Tuple{Int,Int},Tuple{Int,Bool}}() for _ in 1:numfields] - - # `facedict` keeps track of the visited faces. We only need to store the first dof we - # add to the face; if we encounter the same face again we *always* reverse the order. In - # 2D a face (i.e. a line) is uniquely determined by 2 vertices, and in 3D a face (i.e. a - # surface) is uniquely determined by 3 vertices. - facedicts = [Dict{NTuple{dim,Int},Int}() for _ in 1:numfields] - - # Set initial values - nextdof = 1 # next free dof to distribute - - @debug "\n\nCreating dofs\n" - for fh in dh.fieldhandlers - # sort the cellset since we want to loop through the cells in a fixed order - cellnumbers = sort(collect(fh.cellset)) - nextdof = _close!( - dh, - cellnumbers, - dh.field_names, - fh.fields, - nextdof, - vertexdicts, - edgedicts, - facedicts, - ) - end - dh.ndofs[] = maximum(dh.cell_dofs; init=0) - dh.closed[] = true - - return dh, vertexdicts, edgedicts, facedicts - -end - -function _close!(dh::MixedDofHandler{dim}, cellnumbers, global_field_names, fields, nextdof, vertexdicts, edgedicts, facedicts) where {dim} - ip_infos = InterpolationInfo[] - for field in fields - interpolation = field.interpolation - ip_info = InterpolationInfo(interpolation) - push!(ip_infos, ip_info) - # TODO: More than one face dof per face in 3D are not implemented yet. This requires - # keeping track of the relative rotation between the faces, not just the - # direction as for faces (i.e. edges) in 2D. - dim == 3 && @assert !any(x -> x > 1, ip_info.nfacedofs) - end - - # loop over all the cells, and distribute dofs for all the fields - for ci in cellnumbers - @debug "Creating dofs for cell #$ci" - - cell = getcells(getgrid(dh), ci) - len_cell_dofs = length(dh.cell_dofs) - dh.cell_dofs_offset[ci] = len_cell_dofs + 1 - - for (local_num, field) in pairs(fields) - # for (local_num, field_name) in enumerate(field_names) - fi = findfirst(i->i == field.name, global_field_names) - @debug "\tfield: $(field.name)" - ip_info = ip_infos[local_num] - - # Distribute dofs for vertices - nextdof = add_vertex_dofs( - dh.cell_dofs, cell, vertexdicts[fi], field.dim, - ip_info.nvertexdofs, nextdof - ) - - # Distribute dofs for edges (only applicable when dim is 3) - if dim == 3 && (ip_info.dim == 3 || ip_info.dim == 2) - # Regular 3D element or 2D interpolation embedded in 3D space - nentitydofs = ip_info.dim == 3 ? ip_info.nedgedofs : ip_info.nfacedofs - nextdof = add_edge_dofs( - dh.cell_dofs, cell, edgedicts[fi], field.dim, - nentitydofs, nextdof - ) - end - - # Distribute dofs for faces. Filter out 2D interpolations in 3D space, since - # they are added above as edge dofs. - if ip_info.dim == dim - nextdof = add_face_dofs( - dh.cell_dofs, cell, facedicts[fi], field.dim, - ip_info.nfacedofs, nextdof - ) - end - - # Distribute internal dofs for cells - nextdof = add_cell_dofs( - dh.cell_dofs, field.dim, ip_info.ncelldofs, nextdof - ) - end - - # Store number of dofs that were added - dh.ndofs_per_cell[ci] = length(dh.cell_dofs) - len_cell_dofs - - # @debug "Dofs for cell #$ci:\n\t$cell_dofs" - end # cell loop - return nextdof -end - -function add_vertex_dofs(cell_dofs, cell, vertexdict, field_dim, nvertexdofs, nextdof) - for (vi, vertex) in pairs(vertices(cell)) - nvertexdofs[vi] > 0 || continue # skip if no dof on this vertex - @assert nvertexdofs[vi] == 1 - first_dof = vertexdict[vertex] - if first_dof > 0 # reuse dof - for d in 1:field_dim - reuse_dof = first_dof + (d-1) - push!(cell_dofs, reuse_dof) - end - else # create dofs - vertexdict[vertex] = nextdof - for _ in 1:field_dim - push!(cell_dofs, nextdof) - nextdof += 1 - end - end - end - return nextdof -end - -function add_face_dofs(cell_dofs, cell, facedict, field_dim, nfacedofs, nextdof) - @debug @assert all(nfacedofs .<= 1) "Currently only supports interpolations with less that 2 dofs per face" - for (fi, face) in pairs(faces(cell)) - nfacedofs[fi] > 0 || continue # skip if no dof on this face - sface = sortface(face) - token = Base.ht_keyindex2!(facedict, sface) - if token > 0 # haskey(facedict, sface) -> reuse dofs - first_dof = facedict.vals[token] # facedict[sface] - for dof_loc in nfacedofs[fi]:-1:1 # always reverse - for d in 1:field_dim - reuse_dof = first_dof + (d-1) + (dof_loc-1)*field_dim - push!(cell_dofs, reuse_dof) - end - end - else # !haskey(facedict, sface) -> create dofs - Base._setindex!(facedict, nextdof, sface, -token) # facedict[sface] = nextdof - for _ in 1:nfacedofs[fi], _ in 1:field_dim - push!(cell_dofs, nextdof) - nextdof += 1 - end - end - end - return nextdof -end - -function add_edge_dofs(cell_dofs, cell, edgedict, field_dim, nedgedofs, nextdof) - for (ei, edge) in pairs(edges(cell)) - nedgedofs[ei] > 0 || continue # skip if no dof on this edge - sedge, this_dir = sortedge(edge) - token = Base.ht_keyindex2!(edgedict, sedge) - if token > 0 # haskey(edgedict, sedge) -> reuse dofs - first_dof, prev_dir = edgedict.vals[token] # edgedict[sedge] - # For an edge between vertices v1 and v2 with "dof locations" l1 and l2 and - # already distributed dofs d1-d6: v1 --- l1(d1,d2,d3) --- l2(d4,d5,d6) --- v2: - # - this_dir == prev_dir: first_dof is d1 and we loop as usual over dof - # locations then over field dims - # - this_dir != prev_dir: first_dof is d4 and we need to reverse the loop over - # the dof locations but *not* the field dims - for dof_loc in (this_dir == prev_dir ? (1:nedgedofs[ei]) : (nedgedofs[ei]:-1:1)) - for d in 1:field_dim - reuse_dof = first_dof + (d-1) + (dof_loc-1)*field_dim - push!(cell_dofs, reuse_dof) - end - end - else # !haskey(edgedict, sedge) -> create dofs - Base._setindex!(edgedict, (nextdof, this_dir), sedge, -token) # edgedict[sedge] = (nextdof, this_dir) - for _ in 1:nedgedofs[ei], _ in 1:field_dim - push!(cell_dofs, nextdof) - nextdof += 1 - end - end - end - return nextdof -end - -function add_cell_dofs(cell_dofs, field_dim, ncelldofs, nextdof) - for _ in 1:ncelldofs, _ in 1:field_dim - push!(cell_dofs, nextdof) - nextdof += 1 - end - return nextdof -end - -""" - sortedge(edge::Tuple{Int,Int}) - -Returns the unique representation of an edge and its orientation. -Here the unique representation is the sorted node index tuple. The -orientation is `true` if the edge is not flipped, where it is `false` -if the edge is flipped. -""" -function sortedge(edge::Tuple{Int,Int}) - a, b = edge - a < b ? (return (edge, true)) : (return ((b, a), false)) -end - -""" - sortface(face::Tuple{Int,Int}) - sortface(face::Tuple{Int,Int,Int}) - sortface(face::Tuple{Int,Int,Int,Int}) - -Returns the unique representation of a face. -Here the unique representation is the sorted node index tuple. -Note that in 3D we only need indices to uniquely identify a face, -so the unique representation is always a tuple length 3. -""" -sortface(face::Tuple{Int,Int}) = minmax(face[1], face[2]) -function sortface(face::Tuple{Int,Int,Int}) - a, b, c = face - b, c = minmax(b, c) - a, c = minmax(a, c) - a, b = minmax(a, b) - return (a, b, c) -end -function sortface(face::Tuple{Int,Int,Int,Int}) - a, b, c, d = face - c, d = minmax(c, d) - b, d = minmax(b, d) - a, d = minmax(a, d) - b, c = minmax(b, c) - a, c = minmax(a, c) - a, b = minmax(a, b) - return (a, b, c) -end - - -""" - find_field(dh::MixedDofHandler, field_name::Symbol)::NTuple{2,Int} - -Return the index of the field with name `field_name` in a `MixedDofHandler`. The index is a -`NTuple{2,Int}`, where the 1st entry is the index of the `FieldHandler` within which the -field was found and the 2nd entry is the index of the field within the `FieldHandler`. - -!!! note - Always finds the 1st occurence of a field within `MixedDofHandler`. - -See also: [`find_field(fh::FieldHandler, field_name::Symbol)`](@ref), -[`_find_field(fh::FieldHandler, field_name::Symbol)`](@ref). -""" -function find_field(dh::MixedDofHandler, field_name::Symbol) - for (fh_idx, fh) in pairs(dh.fieldhandlers) - field_idx = _find_field(fh, field_name) - !isnothing(field_idx) && return (fh_idx, field_idx) - end - error("Did not find field :$field_name (existing fields: $(getfieldnames(dh))).") -end - -""" - find_field(fh::FieldHandler, field_name::Symbol)::Int - -Return the index of the field with name `field_name` in a `FieldHandler`. Throw an -error if the field is not found. - -See also: [`find_field(dh::MixedDofHandler, field_name::Symbol)`](@ref), [`_find_field(fh::FieldHandler, field_name::Symbol)`](@ref). -""" -function find_field(fh::FieldHandler, field_name::Symbol) - field_idx = _find_field(fh, field_name) - if field_idx === nothing - error("Did not find field :$field_name in FieldHandler (existing fields: $([field.name for field in fh.fields]))") - end - return field_idx -end - -# No error if field not found -""" - _find_field(fh::FieldHandler, field_name::Symbol)::Int - -Return the index of the field with name `field_name` in the `FieldHandler` `fh`. Return -`nothing` if the field is not found. - -See also: [`find_field(dh::MixedDofHandler, field_name::Symbol)`](@ref), [`find_field(fh::FieldHandler, field_name::Symbol)`](@ref). -""" -function _find_field(fh::FieldHandler, field_name::Symbol) - for (field_idx, field) in pairs(fh.fields) - if field.name == field_name - return field_idx - end - end - return nothing -end - -# Calculate the offset to the first local dof of a field -function field_offset(fh::FieldHandler, field_idx::Int) - offset = 0 - for i in 1:(field_idx-1) - offset += getnbasefunctions(fh.fields[i].interpolation)::Int * fh.fields[i].dim - end - return offset -end -field_offset(fh::FieldHandler, field_name::Symbol) = field_offset(fh, find_field(fh, field_name)) - -field_offset(dh::MixedDofHandler, field_name::Symbol) = field_offset(dh, find_field(dh, field_name)) -function field_offset(dh::MixedDofHandler, field_idxs::Tuple{Int, Int}) - fh_idx, field_idx = field_idxs - field_offset(dh.fieldhandlers[fh_idx], field_idx) -end - -""" - dof_range(fh::FieldHandler, field_idx::Int) - dof_range(fh::FieldHandler, field_name::Symbol) - dof_range(dh:MixedDofHandler, field_name::Symbol) - -Return the local dof range for a given field. The field can be specified by its name or -index, where `field_idx` represents the index of a field within a `FieldHandler` and -`field_idxs` is a tuple of the `FieldHandler`-index within the `MixedDofHandler` and the -`field_idx`. - -!!! note - The `dof_range` of a field can vary between different `FieldHandler`s. Therefore, it is - advised to use the `field_idxs` or refer to a given `FieldHandler` directly in case - several `FieldHandler`s exist. Using the `field_name` will always refer to the first - occurence of `field` within the `MixedDofHandler`. - -Example: -```jldoctest -julia> grid = generate_grid(Triangle, (3, 3)) -Grid{2, Triangle, Float64} with 18 Triangle cells and 16 nodes - -julia> dh = MixedDofHandler(grid); add!(dh, :u, 3); add!(dh, :p, 1); close!(dh); - -julia> dof_range(dh, :u) -1:9 - -julia> dof_range(dh, :p) -10:12 - -julia> dof_range(dh, (1,1)) # field :u -1:9 - -julia> dof_range(dh.fieldhandlers[1], 2) # field :p -10:12 -``` -""" -function dof_range(fh::FieldHandler, field_idx::Int) - offset = field_offset(fh, field_idx) - field_interpolation = fh.fields[field_idx].interpolation - field_dim = fh.fields[field_idx].dim - n_field_dofs = getnbasefunctions(field_interpolation)::Int * field_dim - return (offset+1):(offset+n_field_dofs) -end -dof_range(fh::FieldHandler, field_name::Symbol) = dof_range(fh, find_field(fh, field_name)) - -function dof_range(dh::MixedDofHandler, field_name::Symbol) - if length(dh.fieldhandlers) > 1 - error("The given MixedDofHandler has $(length(dh.fieldhandlers)) FieldHandlers. - Extracting the dof range based on the fieldname might not be a unique problem - in this case. Use `dof_range(fh::FieldHandler, field_name)` instead.") - end - fh_idx, field_idx = find_field(dh, field_name) - return dof_range(dh.fieldhandlers[fh_idx], field_idx) -end - -""" - getfieldinterpolation(dh::MixedDofHandler, field_idxs::NTuple{2,Int}) - getfieldinterpolation(dh::FieldHandler, field_idx::Int) - getfieldinterpolation(dh::FieldHandler, field_name::Symbol) - -Return the interpolation of a given field. The field can be specified by its index (see -[`find_field`](@ref) or its name. -""" -function getfieldinterpolation(dh::MixedDofHandler, field_idxs::NTuple{2,Int}) - fh_idx, field_idx = field_idxs - ip = dh.fieldhandlers[fh_idx].fields[field_idx].interpolation - return ip -end -getfieldinterpolation(fh::FieldHandler, field_idx::Int) = fh.fields[field_idx].interpolation -getfieldinterpolation(fh::FieldHandler, field_name::Symbol) = getfieldinterpolation(fh, find_field(fh, field_name)) - -""" - reshape_to_nodes(dh::AbstractDofHandler, u::Vector{T}, fieldname::Symbol) where T - -Reshape the entries of the dof-vector `u` which correspond to the field `fieldname` in nodal order. -Return a matrix with a column for every node and a row for every dimension of the field. -For superparametric fields only the entries corresponding to nodes of the grid will be returned. Do not use this function for subparametric approximations. -""" -function reshape_to_nodes(dh::MixedDofHandler, u::Vector{T}, fieldname::Symbol) where T - # make sure the field exists - fieldname ∈ getfieldnames(dh) || error("Field $fieldname not found.") - - field_dim = getfielddim(dh, fieldname) - space_dim = field_dim == 2 ? 3 : field_dim - data = fill(T(NaN), space_dim, getnnodes(getgrid(dh))) # set default value - - for fh in dh.fieldhandlers - # check if this fh contains this field, otherwise continue to the next - field_idx = _find_field(fh, fieldname) - field_idx === nothing && continue - offset = field_offset(fh, field_idx) - - reshape_field_data!(data, dh, u, offset, field_dim, fh.cellset) - end - return data -end - -function reshape_field_data!(data::Matrix{T}, dh::AbstractDofHandler, u::Vector{T}, field_offset::Int, field_dim::Int, cellset=1:getncells(getgrid(dh))) where T - - for cell in CellIterator(dh, cellset, UpdateFlags(; nodes=true, coords=false, dofs=true)) - _celldofs = celldofs(cell) - counter = 1 - for node in getnodes(cell) - for d in 1:field_dim - data[d, node] = u[_celldofs[counter + field_offset]] - @debug println(" exporting $(u[_celldofs[counter + field_offset]]) for dof#$(_celldofs[counter + field_offset])") - counter += 1 - end - if field_dim == 2 - # paraview requires 3D-data so pad with zero - data[3, node] = 0 - end - end - end - return data -end From 303513e40087e1f0a8219f6d93fd6488e4de76f3 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 10 Jun 2023 11:18:05 +0200 Subject: [PATCH 11/16] getgrid to get_grid --- src/Dofs/ConstraintHandler.jl | 18 +++++++++--------- src/Dofs/DofHandler.jl | 32 ++++++++++++++++---------------- src/Dofs/apply_analytical.jl | 14 +++++++------- src/Export/VTK.jl | 2 +- src/L2_projection.jl | 12 ++++++------ src/exports.jl | 2 +- src/iterators.jl | 12 ++++++------ 7 files changed, 46 insertions(+), 46 deletions(-) diff --git a/src/Dofs/ConstraintHandler.jl b/src/Dofs/ConstraintHandler.jl index 94d6d11fc5..aeb51ed467 100644 --- a/src/Dofs/ConstraintHandler.jl +++ b/src/Dofs/ConstraintHandler.jl @@ -315,8 +315,8 @@ function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, b return local_face_dofs, local_face_dofs_offset end -function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(getgrid(ch.dh)))) - grid = getgrid(ch.dh) +function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(get_grid(ch.dh)))) + grid = get_grid(ch.dh) if interpolation !== default_interpolation(getcelltype(grid, first(cellset))) @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.") end @@ -449,7 +449,7 @@ function _update!(inhomogeneities::Vector{Float64}, f::Function, ::Set{Int}, fie dofmapping::Dict{Int,Int}, dofcoefficients::Vector{Union{Nothing,DofCoefficients{T}}}, time::Real) where T counter = 1 for nodenumber in nodeidxs - x = getcoordinates(getnodes(getgrid(dh), nodenumber)) + x = getcoordinates(getnodes(get_grid(dh), nodenumber)) bc_value = f(x, time) @assert length(bc_value) == length(components) for v in bc_value @@ -477,13 +477,13 @@ function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler) for field in unique_fields nd = getfielddim(ch.dh, field) - data = zeros(Float64, nd, getnnodes(getgrid(ch.dh))) + data = zeros(Float64, nd, getnnodes(get_grid(ch.dh))) for dbc in ch.dbcs dbc.field_name != field && continue if eltype(dbc.faces) <: BoundaryIndex functype = boundaryfunction(eltype(dbc.faces)) for (cellidx, faceidx) in dbc.faces - for facenode in functype(getcells(getgrid(ch.dh), cellidx))[faceidx] + for facenode in functype(getcells(get_grid(ch.dh), cellidx))[faceidx] for component in dbc.components data[component, facenode] = 1 end @@ -823,7 +823,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet) dbc.field_name in fh.field_names || continue # Compute the intersection between dbc.set and the cellset of this # FieldHandler and skip if the set is empty - filtered_set = filter_dbc_set(getgrid(ch.dh), fh.cellset, dbc.faces) + filtered_set = filter_dbc_set(get_grid(ch.dh), fh.cellset, dbc.faces) isempty(filtered_set) && continue # Fetch information about the field on this FieldHandler field_idx = find_field(fh, dbc.field_name) @@ -844,7 +844,7 @@ function add!(ch::ConstraintHandler, dbc::Dirichlet) # BCValues are just dummy for nodesets so set to FaceIndex EntityType = FaceIndex end - CT = getcelltype(getgrid(ch.dh), fh) # Same celltype enforced in FieldHandler constructor + CT = getcelltype(get_grid(ch.dh), fh) # Same celltype enforced in FieldHandler constructor bcvalues = BCValues(interpolation, default_interpolation(CT), EntityType) # Recreate the Dirichlet(...) struct with the filtered set and call internal add! filtered_dbc = Dirichlet(dbc.field_name, filtered_set, dbc.f, components) @@ -940,7 +940,7 @@ function add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet) is_legacy = !isempty(pdbc.face_pairs) && isempty(pdbc.face_map) if is_legacy for (mset, iset) in pdbc.face_pairs - collect_periodic_faces!(pdbc.face_map, getgrid(ch.dh), mset, iset, identity) # TODO: Better transform + collect_periodic_faces!(pdbc.face_map, get_grid(ch.dh), mset, iset, identity) # TODO: Better transform end end field_idx = find_field(ch.dh, pdbc.field_name) @@ -981,7 +981,7 @@ end function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::Interpolation, field_dim::Int, offset::Int, is_legacy::Bool, rotation_matrix::Union{Matrix{T},Nothing}, ::Type{dof_map_t}, iterator_f::F) where {T, dof_map_t, F <: Function} - grid = getgrid(ch.dh) + grid = get_grid(ch.dh) face_map = pdbc.face_map # Indices of the local dofs for the faces diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 3eeadfa99d..900e702276 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -1,7 +1,7 @@ abstract type AbstractDofHandler end """ - getgrid(dh::AbstractDofHandler) + get_grid(dh::AbstractDofHandler) Access the some grid representation for the dof handler. @@ -10,7 +10,7 @@ Access the some grid representation for the dof handler. distributed assembly and assembly on a single process, because most parts oft the functionality can be handled by only acting on the locally owned cell set. """ -getgrid(dh::AbstractDofHandler) +get_grid(dh::AbstractDofHandler) """ Field(name::Symbol, interpolation::Interpolation, dim::Int) @@ -110,7 +110,7 @@ function Base.show(io::IO, ::MIME"text/plain", dh::DofHandler) end isclosed(dh::AbstractDofHandler) = dh.closed[] -getgrid(dh::DofHandler) = dh.grid +get_grid(dh::DofHandler) = dh.grid """ ndofs(dh::AbstractDofHandler) @@ -127,11 +127,11 @@ Return the number of degrees of freedom for the cell with index `cell`. See also [`ndofs`](@ref). """ function ndofs_per_cell(dh::DofHandler, cell::Int=1) - @boundscheck 1 <= cell <= getncells(getgrid(dh)) + @boundscheck 1 <= cell <= getncells(get_grid(dh)) return @inbounds ndofs_per_cell(dh.fieldhandlers[dh.cell_to_fieldhandler[cell]]) end ndofs_per_cell(fh::FieldHandler) = fh.ndofs_per_cell -nnodes_per_cell(dh::DofHandler, cell::Int=1) = nnodes_per_cell(getgrid(dh), cell) # TODO: deprecate, shouldn't belong to DofHandler any longer +nnodes_per_cell(dh::DofHandler, cell::Int=1) = nnodes_per_cell(get_grid(dh), cell) # TODO: deprecate, shouldn't belong to DofHandler any longer """ celldofs!(global_dofs::Vector{Int}, dh::AbstractDofHandler, i::Int) @@ -160,11 +160,11 @@ end #TODO: perspectively remove in favor of `getcoordinates!(global_coords, grid, i)`? function cellcoords!(global_coords::Vector{Vec{dim,T}}, dh::DofHandler, i::Union{Int, <:AbstractCell}) where {dim,T} - cellcoords!(global_coords, getgrid(dh), i) + cellcoords!(global_coords, get_grid(dh), i) end function cellnodes!(global_nodes::Vector{Int}, dh::DofHandler, i::Union{Int, <:AbstractCell}) - cellnodes!(global_nodes, getgrid(dh), i) + cellnodes!(global_nodes, get_grid(dh), i) end """ @@ -204,11 +204,11 @@ Add all fields of the [`FieldHandler`](@ref) `fh` to `dh`. function add!(dh::DofHandler, fh::FieldHandler) # TODO: perhaps check that a field with the same name is the same field? @assert !isclosed(dh) - _check_same_celltype(getgrid(dh), collect(fh.cellset)) + _check_same_celltype(get_grid(dh), collect(fh.cellset)) _check_cellset_intersections(dh, fh) # the field interpolations should have the same refshape as the cells they are applied to # extract the celltype from the first cell as the celltypes are all equal - cell_type = getcelltype(getgrid(dh), fh) + cell_type = getcelltype(get_grid(dh), fh) refshape_cellset = getrefshape(default_interpolation(cell_type)) for interpolation in fh.field_interpolations refshape = getrefshape(interpolation) @@ -235,11 +235,11 @@ celltypes, [`add!(dh::DofHandler, fh::FieldHandler)`](@ref) must be used instead function add!(dh::DofHandler, name::Symbol, ip::Interpolation) @assert !isclosed(dh) - celltype = getcelltype(getgrid(dh)) + celltype = getcelltype(get_grid(dh)) @assert isconcretetype(celltype) if length(dh.fieldhandlers) == 0 - cellset = Set(1:getncells(getgrid(dh))) + cellset = Set(1:getncells(get_grid(dh))) push!(dh.fieldhandlers, FieldHandler(Field[], cellset)) elseif length(dh.fieldhandlers) > 1 error("If you have more than one FieldHandler, you must specify field") @@ -295,7 +295,7 @@ function __close!(dh::DofHandler{dim}) where {dim} # `vertexdict` keeps track of the visited vertices. The first dof added to vertex v is # stored in vertexdict[v]. # TODO: No need to allocate this vector for fields that don't have vertex dofs - vertexdicts = [zeros(Int, getnnodes(getgrid(dh))) for _ in 1:numfields] + vertexdicts = [zeros(Int, getnnodes(get_grid(dh))) for _ in 1:numfields] # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem. # An edge is uniquely determined by two global vertices, with global direction going @@ -393,7 +393,7 @@ function _close_fieldhandler!(dh::DofHandler{sdim}, fh::FieldHandler, fh_index:: @assert dh.cell_to_fieldhandler[ci] == 0 dh.cell_to_fieldhandler[ci] = fh_index - cell = getcells(getgrid(dh), ci) + cell = getcells(get_grid(dh), ci) len_cell_dofs_start = length(dh.cell_dofs) dh.cell_dofs_offset[ci] = len_cell_dofs_start + 1 @@ -835,10 +835,10 @@ function _evaluate_at_grid_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol # VTK output of solution field (or L2 projected scalar data) n_c = n_components(ip) vtk_dim = n_c == 2 ? 3 : n_c # VTK wants vectors padded to 3D - data = fill(NaN * zero(T), vtk_dim, getnnodes(getgrid(dh))) + data = fill(NaN * zero(T), vtk_dim, getnnodes(get_grid(dh))) else # Just evalutation at grid nodes - data = fill(NaN * zero(RT), getnnodes(getgrid(dh))) + data = fill(NaN * zero(RT), getnnodes(get_grid(dh))) end # Loop over the fieldhandlers for fh in dh.fieldhandlers @@ -846,7 +846,7 @@ function _evaluate_at_grid_nodes(dh::DofHandler, u::Vector{T}, fieldname::Symbol field_idx = _find_field(fh, fieldname) field_idx === nothing && continue # Set up CellValues with the local node coords as quadrature points - CT = getcelltype(getgrid(dh), fh) + CT = getcelltype(get_grid(dh), fh) ip_geo = default_interpolation(CT) local_node_coords = reference_coordinates(ip_geo) qr = QuadratureRule{getrefshape(ip)}(zeros(length(local_node_coords)), local_node_coords) diff --git a/src/Dofs/apply_analytical.jl b/src/Dofs/apply_analytical.jl index 6775749228..351b744cf6 100644 --- a/src/Dofs/apply_analytical.jl +++ b/src/Dofs/apply_analytical.jl @@ -1,12 +1,12 @@ function _default_interpolations(dh::DofHandler) fhs = dh.fieldhandlers - ntuple(i -> default_interpolation(getcelltype(getgrid(dh), fhs[i])), length(fhs)) + ntuple(i -> default_interpolation(getcelltype(get_grid(dh), fhs[i])), length(fhs)) end """ apply_analytical!( a::AbstractVector, dh::AbstractDofHandler, fieldname::Symbol, - f::Function, cellset=1:getncells(getgrid(dh))) + f::Function, cellset=1:getncells(get_grid(dh))) Apply a solution `f(x)` by modifying the values in the degree of freedom vector `a` pertaining to the field `fieldname` for all cells in `cellset`. @@ -26,7 +26,7 @@ This function can be used to apply initial conditions for time dependent problem """ function apply_analytical!( a::AbstractVector, dh::DofHandler, fieldname::Symbol, f::Function, - cellset = 1:getncells(getgrid(dh))) + cellset = 1:getncells(get_grid(dh))) fieldname ∉ getfieldnames(dh) && error("The fieldname $fieldname was not found in the dof handler") ip_geos = _default_interpolations(dh) @@ -37,8 +37,8 @@ function apply_analytical!( ip_fun = getfieldinterpolation(fh, field_idx) field_dim = getfielddim(fh, field_idx) celldofinds = dof_range(fh, fieldname) - set_intersection = if length(cellset) == length(fh.cellset) == getncells(getgrid(dh)) - BitSet(1:getncells(getgrid(dh))) + set_intersection = if length(cellset) == length(fh.cellset) == getncells(get_grid(dh)) + BitSet(1:getncells(get_grid(dh))) else intersect(BitSet(fh.cellset), BitSet(cellset)) end @@ -51,7 +51,7 @@ function _apply_analytical!( a::AbstractVector, dh::AbstractDofHandler, celldofinds, field_dim, ip_fun::Interpolation{RefShape}, ip_geo::Interpolation, f::Function, cellset) where {dim, RefShape<:AbstractRefShape{dim}} - coords = getcoordinates(getgrid(dh), first(cellset)) + coords = getcoordinates(get_grid(dh), first(cellset)) ref_points = reference_coordinates(ip_fun) dummy_weights = zeros(length(ref_points)) qr = QuadratureRule{RefShape}(dummy_weights, ref_points) @@ -64,7 +64,7 @@ function _apply_analytical!( length(f(first(coords))) == field_dim || error("length(f(x)) must be equal to dimension of the field ($field_dim)") for cellnr in cellset - getcoordinates!(coords, getgrid(dh), cellnr) + getcoordinates!(coords, get_grid(dh), cellnr) celldofs!(c_dofs, dh, cellnr) for (i, celldofind) in enumerate(celldofinds) f_dofs[i] = c_dofs[celldofind] diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index e88f4f6019..8ddb5e4cef 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -35,7 +35,7 @@ function WriteVTK.vtk_grid(filename::AbstractString, grid::Grid{dim,C,T}; kwargs return vtk_grid(filename, coords, cls; kwargs...) end function WriteVTK.vtk_grid(filename::AbstractString, dh::AbstractDofHandler; kwargs...) - vtk_grid(filename, dh.grid; kwargs...) + vtk_grid(filename, get_grid(dh); kwargs...) end function toparaview!(v, x::Vec{D}) where D diff --git a/src/L2_projection.jl b/src/L2_projection.jl index 95f6dd5bf1..7cfd954c42 100644 --- a/src/L2_projection.jl +++ b/src/L2_projection.jl @@ -73,7 +73,7 @@ _mass_qr(ip::VectorizedInterpolation) = _mass_qr(ip.ip) function Base.show(io::IO, ::MIME"text/plain", proj::L2Projector) println(io, typeof(proj)) - println(io, " projection on: ", length(proj.set), "/", getncells(getgrid(proj.dh)), " cells in grid") + println(io, " projection on: ", length(proj.set), "/", getncells(get_grid(proj.dh)), " cells in grid") println(io, " function interpolation: ", proj.func_ip) println(io, " geometric interpolation: ", proj.geom_ip) end @@ -100,7 +100,7 @@ function _assemble_L2_matrix(fe_values, set, dh) celldofs!(cell_dofs, dh, cellnum) fill!(Me, 0) - Xe = getcoordinates(getgrid(dh), cellnum) + Xe = getcoordinates(get_grid(dh), cellnum) reinit!(fe_values, Xe) ## ∭( v ⋅ u )dΩ @@ -195,7 +195,7 @@ function _project(vars, proj::L2Projector, fe_values::AbstractValues, M::Integer for (ic,cellnum) in enumerate(proj.set) celldofs!(cell_dofs, proj.dh, cellnum) fill!(fe, 0) - Xe = getcoordinates(getgrid(proj.dh), cellnum) + Xe = getcoordinates(get_grid(proj.dh), cellnum) cell_vars = vars[ic] reinit!(fe_values, Xe) @@ -226,7 +226,7 @@ end function WriteVTK.vtk_point_data(vtk::WriteVTK.DatasetFile, proj::L2Projector, vals::Vector{T}, name::AbstractString) where T data = _evaluate_at_grid_nodes(proj, vals, #=vtk=# Val(true))::Matrix - @assert size(data, 2) == getnnodes(getgrid(proj.dh)) + @assert size(data, 2) == getnnodes(get_grid(proj.dh)) vtk_point_data(vtk, data, name; component_names=component_names(T)) return vtk end @@ -248,9 +248,9 @@ function _evaluate_at_grid_nodes( @assert ndofs(dh) == length(vals) if vtk nout = S <: Vec{2} ? 3 : M # Pad 2D Vec to 3D - data = fill(T(NaN), nout, getnnodes(getgrid(dh))) + data = fill(T(NaN), nout, getnnodes(get_grid(dh))) else - data = fill(NaN * zero(S), getnnodes(getgrid(dh))) + data = fill(NaN * zero(S), getnnodes(get_grid(dh))) end ip, gip = proj.func_ip, proj.geom_ip refdim, refshape = getdim(ip), getrefshape(ip) diff --git a/src/exports.jl b/src/exports.jl index 5293d122b5..d61e09183f 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -119,7 +119,7 @@ export Field, evaluate_at_grid_nodes, apply_analytical!, - getgrid, + get_grid, getdim, # Constraints diff --git a/src/iterators.jl b/src/iterators.jl index 39c79c7afc..2040149cf0 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -58,12 +58,12 @@ function CellCache(grid::Grid{dim,C,T}, flags::UpdateFlags=UpdateFlags()) where end function CellCache(dh::DofHandler{dim}, flags::UpdateFlags=UpdateFlags()) where {dim} - N = nnodes_per_cell(getgrid(dh)) + N = nnodes_per_cell(get_grid(dh)) nodes = zeros(Int, N) - coords = zeros(Vec{dim, get_coordinate_eltype(getgrid(dh))}, N) + coords = zeros(Vec{dim, get_coordinate_eltype(get_grid(dh))}, N) n = ndofs_per_cell(dh) celldofs = zeros(Int, n) - return CellCache(flags, getgrid(dh), ScalarWrapper(-1), nodes, coords, dh, celldofs) + return CellCache(flags, get_grid(dh), ScalarWrapper(-1), nodes, coords, dh, celldofs) end function reinit!(cc::CellCache, i::Int) @@ -141,13 +141,13 @@ function CellIterator(gridordh::Union{Grid,AbstractDofHandler}, set::Union{IntegerCollection,Nothing}=nothing, flags::UpdateFlags=UpdateFlags()) if set === nothing - grid = gridordh isa AbstractDofHandler ? getgrid(gridordh) : gridordh + grid = gridordh isa AbstractDofHandler ? get_grid(gridordh) : gridordh set = 1:getncells(grid) end - if gridordh isa DofHandler && !isconcretetype(getcelltype(gridordh.grid)) + if gridordh isa DofHandler && !isconcretetype(getcelltype(getgrid(gridordh))) # TODO: Since the CellCache is resizeable this is not really necessary to check # here, but might be useful to catch slow code paths? - _check_same_celltype(getgrid(gridordh), set) + _check_same_celltype(get_grid(gridordh), set) end return CellIterator(CellCache(gridordh, flags), set) end From 5d663a9a14bfe06b595a0226f3911c4728b7aa60 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 10 Jun 2023 11:23:02 +0200 Subject: [PATCH 12/16] Revert export of getdim. --- src/exports.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/exports.jl b/src/exports.jl index d61e09183f..5412cdbf0b 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -120,7 +120,6 @@ export evaluate_at_grid_nodes, apply_analytical!, get_grid, - getdim, # Constraints ConstraintHandler, From 75a3f97b4b27cc451eba1e03ba77ff3872d48b28 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 10 Jun 2023 11:24:39 +0200 Subject: [PATCH 13/16] Apply suggestions from code review Co-authored-by: Fredrik Ekre --- src/Dofs/DofHandler.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Dofs/DofHandler.jl b/src/Dofs/DofHandler.jl index 900e702276..16945ecd80 100644 --- a/src/Dofs/DofHandler.jl +++ b/src/Dofs/DofHandler.jl @@ -3,11 +3,11 @@ abstract type AbstractDofHandler end """ get_grid(dh::AbstractDofHandler) -Access the some grid representation for the dof handler. +Access some grid representation for the dof handler. !!! note This API function is currently not well-defined. It acts as the interface between - distributed assembly and assembly on a single process, because most parts oft the + distributed assembly and assembly on a single process, because most parts of the functionality can be handled by only acting on the locally owned cell set. """ get_grid(dh::AbstractDofHandler) From 0f4908de0c8028624d0d37a59f05bb7243643c97 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 10 Jun 2023 11:29:46 +0200 Subject: [PATCH 14/16] Docs --- docs/src/devdocs/dofhandler.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/devdocs/dofhandler.md b/docs/src/devdocs/dofhandler.md index b4dff58cf0..531def5284 100644 --- a/docs/src/devdocs/dofhandler.md +++ b/docs/src/devdocs/dofhandler.md @@ -19,7 +19,7 @@ Ferrite.SurfaceOrientationInfo The main entry point for dof distribution is [`__close!`](@ref). ```@docs -Ferrite.getgrid +Ferrite.get_grid Ferrite.find_field(dh::DofHandler, field_name::Symbol) Ferrite._find_field(fh::FieldHandler, field_name::Symbol) Ferrite._close_fieldhandler! From 5951b79ab4a3edf2e2336df4f2615f3810f34ff7 Mon Sep 17 00:00:00 2001 From: Dennis Ogiermann Date: Sat, 10 Jun 2023 11:30:44 +0200 Subject: [PATCH 15/16] Oopsie. --- src/iterators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/iterators.jl b/src/iterators.jl index 2040149cf0..ccc0c8933b 100644 --- a/src/iterators.jl +++ b/src/iterators.jl @@ -144,7 +144,7 @@ function CellIterator(gridordh::Union{Grid,AbstractDofHandler}, grid = gridordh isa AbstractDofHandler ? get_grid(gridordh) : gridordh set = 1:getncells(grid) end - if gridordh isa DofHandler && !isconcretetype(getcelltype(getgrid(gridordh))) + if gridordh isa DofHandler && !isconcretetype(getcelltype(get_grid(gridordh))) # TODO: Since the CellCache is resizeable this is not really necessary to check # here, but might be useful to catch slow code paths? _check_same_celltype(get_grid(gridordh), set) From b65284ff6337d0fc7ce86c5dd3cb596d16fb7893 Mon Sep 17 00:00:00 2001 From: termi-official Date: Tue, 13 Jun 2023 13:54:30 +0200 Subject: [PATCH 16/16] Remote get_grid export. --- src/exports.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/exports.jl b/src/exports.jl index 5412cdbf0b..627b3f32c6 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -119,7 +119,6 @@ export Field, evaluate_at_grid_nodes, apply_analytical!, - get_grid, # Constraints ConstraintHandler,