Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

bugfixes, ntriangles for higher order triangle and introduction of new kwarg copy_fields in interpolategradient() #83

Merged
merged 9 commits into from
May 24, 2023
21 changes: 12 additions & 9 deletions src/makieplotting.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ function Makie.plot!(SP::SolutionPlot{<:Tuple{<:MakiePlotter}})
mins = @lift(minimum(x->isnan(x) ? 1e10 : x, $solution))
maxs = @lift(maximum(x->isnan(x) ? -1e10 : x, $solution))
SP[:colorrange] = @lift(isapprox($mins,$maxs) ? ($mins,1.01($maxs)) : ($mins,$maxs))
return Makie.mesh!(SP, plotter.mesh, color=solution, shading=SP[:shading], scale_plot=SP[:scale_plot], colormap=SP[:colormap], transparent=SP[:transparent])
return Makie.mesh!(SP, plotter.mesh, color=solution, shading=SP[:shading], scale_plot=SP[:scale_plot], colormap=SP[:colormap],colorrange=SP[:colorrange] , transparent=SP[:transparent])
end

"""
Expand Down Expand Up @@ -174,9 +174,10 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
# u_matrix = @lift($(WF[:deformation_field])===:default ? zeros(0,3) : transfer_solution(plotter; field_idx=Ferrite.find_field(plotter.dh,$(WF[:deformation_field])), process=identity))
# coords = @lift($(WF[:deformation_field])===:default ? plotter.physical_coords : plotter.physical_coords .+ ($(WF[:scale]) .* $(u_matrix)))
#original representation
nodal_u_matrix = @lift begin
pointtype = dim > 2 ? Point3f : Point2f
termi-official marked this conversation as resolved.
Show resolved Hide resolved
nodal_u_matrix = @lift begin
if $(WF[:deformation_field])===:default
Point3f[Point3f(0,0,0)]
pointtype[zero(pointtype)]
else
convert(Vector{Point{Ferrite.getdim(plotter.dh.grid),Float32}},Makie.to_vertices(dof_to_node(plotter.dh, $(WF[1][].u); field_name=$(WF[:deformation_field]))))
end
Expand All @@ -196,7 +197,7 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
end
lines
end
nodes = @lift($(WF[:plotnodes]) ? $(gridnodes) : Point3f[Point3f(0,0,0)])
nodes = @lift($(WF[:plotnodes]) ? $(gridnodes) : pointtype[zero(pointtype)])
#plot cellsets
cellsets = plotter.dh.grid.cellsets
cellset_to_value = Dict{String,Int}()
Expand All @@ -213,9 +214,9 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
end
u_matrix = @lift begin
if $(WF[:deformation_field])===:default
Point3f[Point3f(0,0,0)]
pointtype[zero(pointtype)]
else
convert(Vector{Point{Ferrite.getdim(plotter.dh.grid),Float32}},Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_name=$(WF[:deformation_field]), process=identity)))
Makie.to_ndim.(pointtype, Makie.to_vertices(transfer_solution(plotter,$(plotter.u); field_idx=Ferrite.find_field(plotter.dh,$(WF[:deformation_field])), process=identity)), 0f0)
end
end
coords = @lift begin
Expand All @@ -229,7 +230,8 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:MakiePlotter{dim}}}) where dim
cellset_u = reshape(transfer_scalar_celldata(plotter, cellset_u; process=identity), num_vertices(plotter))
Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
#plot the nodes
Makie.scatter!(WF,gridnodes,markersize=WF[:markersize], color=WF[:color], visible=WF[:visible])
shouldplot = @lift ($(WF[:visible]) && $(WF[:plotnodes]))
Makie.scatter!(WF,gridnodes,markersize=WF[:markersize], color=WF[:color], visible=shouldplot)
#set up nodelabels
nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size($gridnodes,1)] : [""]
nodepositions = @lift $(WF[:nodelabels]) ? $gridnodes : (dim < 3 ? Point2f[Point2f((0,0))] : Point3f[Point3f((0,0,0))])
Expand All @@ -253,7 +255,8 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where
append!(lines, [coords[e] for boundary in boundaryentities for e in boundary])
end
nodes = @lift($(WF[:plotnodes]) ? coords : Point3f[Point3f(0,0,0)])
Makie.scatter!(WF,nodes,markersize=WF[:markersize], color=WF[:color])
shouldplot = @lift ($(WF[:visible]) && $(WF[:plotnodes]))
Makie.scatter!(WF,nodes,markersize=WF[:markersize], color=WF[:color], visible=shouldplot)
nodelabels = @lift $(WF[:nodelabels]) ? ["$i" for i in 1:size(coords,1)] : [""]
nodepositions = @lift $(WF[:nodelabels]) ? coords : (dim < 3 ? Point2f[Point2f((0,0))] : Point3f[Point3f((0,0,0))])
celllabels = @lift $(WF[:celllabels]) ? ["$i" for i in 1:Ferrite.getncells(grid)] : [""]
Expand Down Expand Up @@ -283,7 +286,7 @@ function Makie.plot!(WF::Wireframe{<:Tuple{<:Ferrite.AbstractGrid{dim}}}) where
Makie.mesh!(WF, plotter.mesh, color=cellset_u, shading=false, scale_plot=false, colormap=:darktest, visible=WF[:cellsets])
Makie.text!(WF,nodelabels, position=nodepositions, textsize=WF[:textsize], offset=WF[:offset],color=WF[:nodelabelcolor])
Makie.text!(WF,celllabels, position=cellpositions, textsize=WF[:textsize], color=WF[:celllabelcolor], align=(:center,:center))
Makie.linesegments!(WF,lines,color=WF[:color], strokewidth=WF[:strokewidth])
Makie.linesegments!(WF,lines,color=WF[:color], strokewidth=WF[:strokewidth], visible=WF[:visible])
end

"""
Expand Down
23 changes: 18 additions & 5 deletions src/utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,7 @@ num_vertices(p::MakiePlotter) = size(p.physical_coords,1)
TODO this looks faulty...think harder.
"""
# Helper to count triangles e.g. for preallocations.
ntriangles(cell::Ferrite.AbstractCell{2,3,3}) = 1 # Tris in 2D
ntriangles(cell::Ferrite.AbstractCell{2,N,3}) where {N} = 1 # Tris in 2D
termi-official marked this conversation as resolved.
Show resolved Hide resolved
ntriangles(cell::Ferrite.AbstractCell{3,3,1}) = 1 # Tris in 3D
ntriangles(cell::Ferrite.AbstractCell{dim,N,4}) where {dim,N} = 4 # Quads in 2D and 3D
ntriangles(cell::Ferrite.AbstractCell{3,N,1}) where N = 4 # Tets as a special case of a Quad, obviously :)
Expand Down Expand Up @@ -452,11 +452,13 @@ This is a helper to access the correct value in Tensors.jl entities, because the
@inline _tensorsjl_gradient_accessor(m::Tensors.Tensor{2,dim}, field_dim_idx::Int, spatial_dim_idx::Int) where {dim} = m[field_dim_idx, spatial_dim_idx]

"""
interpolate_gradient_field(dh::DofHandler, u::AbstractVector, field_name::Symbol)
interpolate_gradient_field(dh::DofHandler, u::AbstractVector, field_name::Symbol; copy_fields::Vector{Symbol})

Compute the piecewise discontinuous gradient field for `field_name`. Returns the flux dof handler and the corresponding flux dof values.
If the additional keyword argument `copy_fields` is provided with a non empty `Vector{Symbol}`, the corresponding fields of `dh` will be
copied into the returned flux dof handler and flux dof value vector.
"""
function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::AbstractVector, field_name::Symbol) where {spatial_dim}
function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::AbstractVector, field_name::Symbol; copy_fields::Vector{Symbol}=Symbol[]) where {spatial_dim}
# Get some helpers
field_idx = Ferrite.find_field(dh, field_name)
ip = Ferrite.getfieldinterpolation(dh, field_idx)
Expand All @@ -466,6 +468,12 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
ip_gradient = get_gradient_interpolation(ip)
field_dim = Ferrite.getfielddim(dh,field_name)
push!(dh_gradient, :gradient, field_dim*spatial_dim, ip_gradient) # field dim × spatial dim components
for fieldname in copy_fields
_field_idx = Ferrite.find_field(dh, fieldname)
_ip = Ferrite.getfieldinterpolation(dh, _field_idx)
_field_dim = Ferrite.getfielddim(dh,fieldname)
push!(dh_gradient, fieldname, _field_dim, _ip)
end
Ferrite.close!(dh_gradient)

num_base_funs = Ferrite.getnbasefunctions(ip_gradient)
Expand All @@ -483,7 +491,7 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst
# Allocate storage for the fluxes to store
u_gradient = zeros(Ferrite.ndofs(dh_gradient))
# In general uᵉ_gradient is an order 3 tensor [field_dim, spatial_dim, num_base_funs]
uᵉ_gradient = zeros(length(cell_dofs_gradient))
uᵉ_gradient = zeros(length(cell_dofs_gradient[Ferrite.dof_range(dh_gradient, :gradient)]))
uᵉ_gradient_view = reshape(uᵉ_gradient, (spatial_dim, field_dim, num_base_funs))
uᵉ = zeros(field_dim*Ferrite.getnbasefunctions(ip))

Expand All @@ -508,7 +516,12 @@ function interpolate_gradient_field(dh::Ferrite.DofHandler{spatial_dim}, u::Abst

# We finally write back the result to the global dof vector of the gradient field
Ferrite.celldofs!(cell_dofs_gradient, dh_gradient, cell_num)
u_gradient[cell_dofs_gradient] .+= uᵉ_gradient
u_gradient[cell_dofs_gradient[Ferrite.dof_range(dh_gradient, :gradient)]] .+= uᵉ_gradient

# copy all requested fields
for fieldname in copy_fields
u_gradient[cell_dofs_gradient[Ferrite.dof_range(dh_gradient, fieldname)]] .= u[cell_dofs[Ferrite.dof_range(dh, fieldname)]]
end
end
return dh_gradient, u_gradient
end
Expand Down