Skip to content

Commit

Permalink
Prevent ambiguity
Browse files Browse the repository at this point in the history
  • Loading branch information
j-fu committed Oct 22, 2024
1 parent bc6dee9 commit 930e76b
Showing 1 changed file with 30 additions and 33 deletions.
63 changes: 30 additions & 33 deletions ext/VoronoiFVMExtendableFEMBaseExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ $(SIGNATURES)
Compute [`VoronoiFVM.edgevelocities`](@ref) for a finite element flow field computed
by [`ExtendableFEM`](https://github.com/chmerdon/ExtendableFEM.jl).
"""
function VoronoiFVM.edgevelocities(grid, vel::FEVectorBlock; kwargs...)
function VoronoiFVM.edgevelocities(grid::ExtendableGrid, vel::FEVectorBlock; kwargs...)
# construct an augmented type to gather precomputed information
# in order to pass it to the repeated integrate call in VoronoiFVM.edgevelocities
axisymmetric = grid[CoordinateSystem] <: Cylindrical2D ? true : false
Expand All @@ -133,7 +133,7 @@ $(SIGNATURES)
Compute [`VoronoiFVM.bfacevelocities`](@ref) for a finite element flow field computed
by [`ExtendableFEM`](https://github.com/chmerdon/ExtendableFEM.jl).
"""
function VoronoiFVM.bfacevelocities(grid, vel::FEVectorBlock; kwargs...)
function VoronoiFVM.bfacevelocities(grid::ExtendableGrid, vel::FEVectorBlock; kwargs...)
axisymmetric = grid[CoordinateSystem] <: Cylindrical2D ? true : false
seg_integrator, point_evaluator, cf, flowgrid = prepare_segment_integration(vel; axisymmetric, kwargs...)
aug_fevec_block = AugmentedFEVectorBlock(vel, seg_integrator, point_evaluator, cf, flowgrid)
Expand Down Expand Up @@ -161,19 +161,16 @@ end
# compute the path integral for the velocity in aug_vec_block between p1 and p2 by
# incrementally walking through each cell in the grid between p1 and p2
# and summing up each cell's contribution
function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVectorBlock; interpolate_eps=1.0e-12, axisymmetric=false, kwargs...)


function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVectorBlock{TVB, TSE, TPE, TCF, TFG}; interpolate_eps=1.0e-12, axisymmetric=false, kwargs...) where {TVB, TSE, TPE, TCF, TFG}
edge_length = dist(p1,p2)
avg_r = (p1[1] + p2[1]) / 2
(; bp1, result, summand, bp2, pint, bpint, bary)= aug_vec_block


if axisymmetric && avg_r < eps()
return 0
end

CF = aug_vec_block.cellfinder
icell::Int = gFindLocal!(bp1, CF, p1; eps=interpolate_eps)
if edge_length interpolate_eps
Expand All @@ -185,10 +182,10 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe
return dot(p2, hnormal)
end
end

SI = aug_vec_block.seg_integrator


xCellFaces::Matrix{Int} = CF.xCellFaces
xFaceCells::Matrix{Int} = CF.xFaceCells
facetogo = CF.facetogo
Expand All @@ -198,41 +195,41 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe
L2G = CF.L2G4EG[1]
L2Gb::Vector{Float64} = L2G.b
invA = CF.invA

t = 0.0

p1_temp2 = @MVector zeros(2)
p1_temp3 = @MVector zeros(3)

function calc_barycentric_coords!(bp, p)
for j = 1:2
cx[j] = p[j] - L2Gb[j]
cx[j] = p[j] - L2Gb[j]
end

fill!(bp, 0)
for j = 1:2, k = 1:2
bp[k] += invA[j, k] * cx[j]
end
postprocess_xreftest!(bp, CF.xCellGeometries[icell])
end

while (true)

# TODO implement proper emergency guard to avoid indefinite loops

# first compute the barycentric coordinates of
# p1,p2

# update local 2 global map
L2G = CF.L2G4EG[1]
update_trafo!(L2G, icell)
L2Gb = L2G.b

mapderiv!(invA, L2G, p1)
calc_barycentric_coords!(bp1, p1)
calc_barycentric_coords!(bp2, p2)

mapderiv!(invA, L2G, p1)

calc_barycentric_coords!(bp1, p1)
calc_barycentric_coords!(bp2, p2)
# if p1 is a node of a triangle, start with
# a cell containing p1 in the direction of (p2-p1)
if count(<=(interpolate_eps), bp1) == 2 # 10^(-13)
Expand All @@ -248,20 +245,20 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe
continue
end
end

# push p1 a little towards the triangle circumcenter
# to avoid it being situated on an edge or node of the
# flowgrid
# (to avoid being stuck on an edge if (p2-p1) is
# on the edge
@. bp1 += interpolate_eps * (bary - bp1)
@. bp1 += interpolate_eps * (bary - bp1)
eval_trafo!(p1, L2G, bp1)

(λp2min, imin) = findmin(bp2)

# if λp2min≥0, then p2 is inside icell and we can simply add the
# integral across the line segment between p1 and p2 to the result

# if not, then p2 is outside of icell and we try to determine
# pint which is the point where icell intersects the line segment
# [p1,p2] - since pint should be on the boundary of icell,
Expand All @@ -270,7 +267,7 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe
# this is not necessarily the previous imin and we have to check
# all triangle edges for if going towards that edge actually takes us
# closer to p2

if λp2min -interpolate_eps
SI.integrator(summand, ((p1, p2)), (bp1, bp2), icell) #!!! allocates; probably due to .integrator being Any
result += summand
Expand Down Expand Up @@ -306,12 +303,12 @@ function _integrate_along_segments(p1, p2, hnormal, aug_vec_block::AugmentedFEVe
eval_trafo!(pint, L2G, bpint)
SI.integrator(summand, ((p1, pint)), (bp1, bpint), icell)
result += summand

# proceed to next cell along edge of smallest barycentric coord
prevcell = icell
icell = xFaceCells[1, xCellFaces[facetogo[1][imin], icell]] #!!! allocates
icell = icell == prevcell ? xFaceCells[2, xCellFaces[facetogo[1][imin], icell]] : icell #!!! allocates

p1 .= pint
end
end
Expand Down

0 comments on commit 930e76b

Please sign in to comment.