Refactor to use grid and dofhandler interface internally

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,21 +1098,22 @@ 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]))
         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)
             cellidx, faceidx = faceidx
             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)
 
@@ -1386,20 +1387,20 @@ 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 = get_coordinate_type(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
     image_mean_x = Tx[]
     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
@@ -1476,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(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
@@ -1511,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, j)) for j 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
@@ -1530,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[i]))
         xij = xj - xi
         if norm(xij - xmij) < TOLh
             found = true
@@ -1546,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
@@ -1593,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

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)
 

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)
 
@@ -645,7 +652,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)

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]

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

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="")