Skip to content
Snippets Groups Projects
Commit b7086df1 authored by Stephan Hilb's avatar Stephan Hilb
Browse files

improve sampling

parent 4d8718da
No related branches found
No related tags found
No related merge requests found
using Statistics: mean
using StaticArrays: SA, SArray, MVector, MMatrix
using StaticArrays: SA, SArray, MVector, MMatrix, SUnitRange
using DataFrames, CSV
export FeSpace, Mapper, FeFunction, ImageFunction, P1, DP0, DP1
export interpolate!, sample, bind!, evaluate, integrate, nabla, save_csv
......@@ -235,6 +235,7 @@ end
Base.show(io::IO, ::MIME"text/plain", x::Derivative) =
print("$(nameof(typeof(x))) of $(typeof(x.f))")
# TODO: should be called "derivative"
function nabla(f)
d = ndims_domain(f.space.mesh)
Derivative(f, zero(MMatrix{d, d, Float64}))
......@@ -317,17 +318,22 @@ end
bind!(f::FacetDivergence, cell) = f.cell[] = cell
function sample(f::FeFunction)
out = _sample(f)
return f.space.size == (1,) ?
reshape(out, size(out)[2:end]) : out
end
function _sample(f::FeFunction)
# assumes dual grid
mesh = f.space.mesh
m0 = NTuple{2}(minimum(mesh.vertices, dims = 2))
m1 = NTuple{2}(maximum(mesh.vertices, dims = 2))
# TODO: make typesafe when introducing zero-dim space size ()
outsize = f.space.size == (1,) ?
round.(Int, m1 .- m0 .+ 1) : (f.space.size, round.(Int, m1 .- m0 .+ 1)...)
fcolon = ntuple(_ -> :, length(f.space.size))
outsize = (f.space.size..., round.(Int, m1 .- m0 .+ 1)...)
out = similar(f.data, outsize)
for cell in cells(mesh)
bind!(f, cell)
......@@ -341,9 +347,11 @@ function sample(f::FeFunction)
for I in CartesianIndex(I0[1], I0[2]):CartesianIndex(I1[1], I1[2])
x = SVector(Tuple(I) .+ m0)
xloc = (A[:, 2:3] .- A[:, 1]) \ (x .- A[:, 1])
# TODO: move to global-to-local function or inside-triangle test
xloc = (A[:, SUnitRange(2, 3)] .- A[:, 1])::SArray \
(x .- A[:, 1])
if all(xloc .>= 0) && sum(xloc) <= 1
out[I + CartesianIndex(1, 1)] = evaluate(f, xloc)[1]
out[fcolon..., I + CartesianIndex(1, 1)] .= evaluate(f, xloc)
end
end
end
......
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment