Merge master

docs/src/devdocs/interpolations.md

@@ -15,6 +15,9 @@ Ferrite.shape_gradient(::Interpolation, ::Vec, ::Int)
 Ferrite.shape_gradient_and_value
 Ferrite.boundarydof_indices
 Ferrite.dirichlet_boundarydof_indices
+Ferrite.shape_values!
+Ferrite.shape_gradients!
+Ferrite.shape_gradients_and_values!
 ```
 
 ### Required methods to implement for all subtypes of `Interpolation` to define a new finite element

docs/src/reference/dofhandler.md

@@ -12,6 +12,7 @@ SubDofHandler
 ## Adding fields to the DofHandlers
 ```@docs
 add!(::DofHandler, ::Symbol, ::Interpolation)
+add!(::SubDofHandler, ::Symbol, ::Interpolation)
 close!(::DofHandler)
 ```
 

src/Dofs/DofHandler.jl

@@ -12,7 +12,6 @@ Access some grid representation for the dof handler.
 """
 get_grid(dh::AbstractDofHandler)
 
-
 struct SubDofHandler{DH} <: AbstractDofHandler
     # From constructor
     dh::DH
@@ -26,8 +25,36 @@ struct SubDofHandler{DH} <: AbstractDofHandler
     # const dof_ranges::Vector{UnitRange{Int}} # TODO: Why not?
 end
 
-# TODO: Should be an inner constructor.
+"""
+    SubDofHandler(dh::AbstractDofHandler, cellset::Set{Int})
+
+Create an `sdh::SubDofHandler` from the parent `dh`, pertaining to the 
+cells in `cellset`. This allows you to add fields to parts of the domain, or using 
+different interpolations or cell types (e.g. `Triangles` and `Quadrilaterals`). All 
+fields and cell types must be the same in one `SubDofHandler`.
+
+After construction any number of discrete fields can be added to the SubDofHandler using
+[`add!`](@ref). Construction is finalized by calling [`close!`](@ref) on the parent `dh`.
+
+# Examples
+We assume we have a `grid` containing "Triangle" and "Quadrilateral" cells, 
+including the cellsets "triangles" and "quadilaterals" for to these cells. 
+```julia
+dh = DofHandler(grid)
+
+sdh_tri = SubDofHandler(dh, getcellset(grid, "triangles"))
+ip_tri = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u
+add!(sdh_tri, :u, ip_tri)
+
+sdh_quad = SubDofHandler(dh, getcellset(grid, "quadilaterals"))
+ip_quad = Lagrange{RefQuadrilateral, 2}()^2 # vector interpolation for a field u
+add!(sdh_quad, :u, ip_quad)
+
+close!(dh) # Finalize by closing the parent 
+```
+"""
 function SubDofHandler(dh::DH, cellset) where {DH <: AbstractDofHandler}
+    # TODO: Should be an inner constructor.
     isclosed(dh) && error("DofHandler already closed")
     # Compute the celltype and make sure all elements have the same one
     CT = getcelltype(dh.grid, first(cellset))
@@ -66,13 +93,6 @@ function _print_field_information(io::IO, mime::MIME"text/plain", sdh::SubDofHan
     end
 end
 
-"""
-    DofHandler(grid::Grid)
-
-Construct a `DofHandler` 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 DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     subdofhandlers::Vector{SubDofHandler{DofHandler{dim, G}}}
     field_names::Vector{Symbol}
@@ -86,6 +106,28 @@ struct DofHandler{dim,G<:AbstractGrid{dim}} <: AbstractDofHandler
     ndofs::ScalarWrapper{Int}
 end
 
+"""
+    DofHandler(grid::Grid)
+
+Construct a `DofHandler` based on the grid `grid`.
+
+After construction any number of discrete fields can be added to the DofHandler using
+[`add!`](@ref). Construction is finalized by calling [`close!`](@ref).
+
+By default fields are added to all elements of the grid. Refer to [`SubDofHandler`](@ref)
+for restricting fields to subdomains of the grid.
+
+# Examples
+
+```julia
+dh = DofHandler(grid)
+ip_u = Lagrange{RefTriangle, 2}()^2 # vector interpolation for a field u
+ip_p = Lagrange{RefTriangle, 1}()   # scalar interpolation for a field p
+add!(dh, :u, ip_u)
+add!(dh, :p, ip_p)
+close!(dh)
+```
+"""
 function DofHandler(grid::G) where {dim, G <: AbstractGrid{dim}}
     ncells = getncells(grid)
     sdhs = SubDofHandler{DofHandler{dim, G}}[]

src/FEValues/FunctionValues.jl

@@ -71,11 +71,8 @@ function Base.copy(v::FunctionValues)
 end
 
 function precompute_values!(fv::FunctionValues, qr::QuadratureRule)
-    n_shape = getnbasefunctions(fv.ip)
     for (qp, ξ) in pairs(getpoints(qr))
-        for i in 1:n_shape
-            fv.dNdξ[i, qp], fv.N_ξ[i, qp] = shape_gradient_and_value(fv.ip, ξ, i)
-        end
+        shape_gradients_and_values!(@view(fv.dNdξ[:, qp]), @view(fv.N[:, qp]), fv.ip, ξ)
     end
 end
 

src/FEValues/GeometryMapping.jl

@@ -50,28 +50,38 @@ function GeometryMapping(::Type{T}, ip::ScalarInterpolation, qr::QuadratureRule,
     VT   = Vec{getdim(ip),T}
     M    = zeros(T,  n_shape, n_qpoints)
     dMdξ = zeros(VT, n_shape, n_qpoints)
-    if RH
-        HT = Tensor{2,getdim(ip),T}
-        dM2dξ2 = zeros(HT, n_shape, n_qpoints)
-    else
-        dM2dξ2 = nothing
-    end
-    for (qp, ξ) in pairs(getpoints(qr))
-        for i in 1:n_shape
-            if RH
-                dM2dξ2[i, qp], dMdξ[i, qp], M[i, qp] = shape_hessian_gradient_and_value(ip, ξ, i)
-            else
-                dMdξ[i, qp], M[i, qp] = shape_gradient_and_value(ip, ξ, i)
-            end
-        end
-    end
+
+    # Hessian (if needed, else nothing)
+    HT = Tensor{2,getdim(ip),T}
+    dM2dξ2 = RH ? zeros(HT, n_shape, n_qpoints) : nothing
+
+    geo_mapping = GeometryMapping(ip, M, dMdξ, dM2dξ2)
+    precompute_values!(geo_mapping, qr)
     return GeometryMapping(ip, M, dMdξ, dM2dξ2)
 end
 function Base.copy(v::GeometryMapping)
     d2Mdξ2_copy = v.d2Mdξ2 === nothing ? nothing : copy(v.d2Mdξ2)
     return GeometryMapping(copy(v.ip), copy(v.M), copy(v.dMdξ), d2Mdξ2_copy)
 end
 
+precompute_values!(gm::GeometryMapping, qr::QuadratureRule) = precompute_values!(gm, RequiresHessian(gm), qr)
+
+function precompute_values!(gm::GeometryMapping, ::RequiresHessian{false}, qr::QuadratureRule)
+    ip = get_geometric_interpolation(gm)
+    for (qp, ξ) in pairs(getpoints(qr))
+        shape_gradients_and_values!(@view(gm.dMdξ[:, qp]), @view(gm.M[:, qp]), ip, ξ)
+    end
+end
+
+function precompute_values!(gm::GeometryMapping, ::RequiresHessian{true}, qr::QuadratureRule)
+    ip = get_geometric_interpolation(gm)
+    for (qp, ξ) in pairs(getpoints(qr))
+        for i in 1:getngeobasefunctions(gm)
+            gm.dM2dξ2[i, qp], gm.dMdξ[i, qp], gm.M[i, qp] = shape_hessian_gradient_and_value(ip, ξ, i)
+        end
+    end
+end
+
 getngeobasefunctions(geo_mapping::GeometryMapping) = size(geo_mapping.M, 1)
 @propagate_inbounds geometric_value(geo_mapping::GeometryMapping, q_point::Int, base_func::Int) = geo_mapping.M[base_func, q_point]
 get_geometric_interpolation(geo_mapping::GeometryMapping) = geo_mapping.ip

src/PointEval/PointEvalHandler.jl

@@ -125,6 +125,7 @@ function find_local_coordinate(interpolation, cell_coordinates::Vector{V}, globa
     for _ in 1:max_iters
         global_guess = zero(V)
         J = zero(Tensor{2, dim, T})
+        # TODO batched eval after 764 is merged.
         for j in 1:n_basefuncs
             dNdξ, N = shape_gradient_and_value(interpolation, local_guess, j)
             global_guess += N * cell_coordinates[j]

src/PointEval/point_values.jl

@@ -91,10 +91,7 @@ shape_value(pv::PointValuesInternal, qp::Int, i::Int) = (@assert qp == 1; pv.N[i
 
 # allow on-the-fly updating
 function reinit!(pv::PointValuesInternal{IP}, coord::Vec{dim}) where {dim, shape <: AbstractRefShape{dim}, IP <: Interpolation{shape}}
-    n_func_basefuncs = getnbasefunctions(pv.ip)
-    for i in 1:n_func_basefuncs
-        pv.N[i] = shape_value(pv.ip, coord, i)
-    end
+    shape_values!(pv.N, pv.ip, coord)
     return nothing
 end
 

src/interpolations.jl

@@ -203,6 +203,46 @@ getnbasefunctions(::Interpolation)
 #   edgedof: dof on a line between 2 vertices (i.e. "corners") (3D only)
 #   celldof: dof that is local to the element
 
+"""
+    shape_values!(values::AbstractArray{T}, ip::Interpolation, ξ::Vec)
+
+Evaluate all shape functions of `ip` at once at the reference point `ξ` and store them in
+`values`.
+"""
+@propagate_inbounds function shape_values!(values::AT, ip::IP, ξ::Vec) where {IP <: Interpolation, AT <: AbstractArray}
+    @boundscheck checkbounds(values, 1:getnbasefunctions(ip))
+    @inbounds for i in 1:getnbasefunctions(ip)
+        values[i] = shape_value(ip, ξ, i)
+    end
+end
+
+"""
+    shape_gradients!(gradients::AbstractArray, ip::Interpolation, ξ::Vec)
+
+Evaluate all shape function gradients of `ip` at once at the reference point `ξ` and store
+them in `gradients`.
+"""
+function shape_gradients!(gradients::AT, ip::IP, ξ::Vec) where {IP <: Interpolation, AT <: AbstractArray}
+    @boundscheck checkbounds(gradients, 1:getnbasefunctions(ip))
+    @inbounds for i in 1:getnbasefunctions(ip)
+        gradients[i] = shape_gradient(ip, ξ, i)
+    end
+end
+
+"""
+    shape_gradients_and_values!(gradients::AbstractArray, values::AbstractArray, ip::Interpolation, ξ::Vec)
+
+Evaluate all shape function gradients and values of `ip` at once at the reference point `ξ`
+and store them in `values`.
+"""
+function shape_gradients_and_values!(gradients::GAT, values::SAT, ip::IP, ξ::Vec) where {IP <: Interpolation, SAT <: AbstractArray, GAT <: AbstractArray}
+    @boundscheck checkbounds(gradients, 1:getnbasefunctions(ip))
+    @boundscheck checkbounds(values, 1:getnbasefunctions(ip))
+    @inbounds for i in 1:getnbasefunctions(ip)
+        gradients[i], values[i] = shape_gradient_and_value(ip, ξ, i)
+    end
+end
+
 """
     shape_value(ip::Interpolation, ξ::Vec, i::Int)
 

test/test_cellvalues.jl

@@ -1,5 +1,5 @@
 @testset "CellValues" begin
-for (scalar_interpol, quad_rule) in  (
+@testset "ip=$scalar_interpol quad_rule=$(typeof(quad_rule))" for (scalar_interpol, quad_rule) in  (
                                     (Lagrange{RefLine, 1}(), QuadratureRule{RefLine}(2)),
                                     (Lagrange{RefLine, 2}(), QuadratureRule{RefLine}(2)),
                                     (Lagrange{RefQuadrilateral, 1}(), QuadratureRule{RefQuadrilateral}(2)),