From 930e76bddc2426811c2f12914793b94993b238e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 22 Oct 2024 16:15:11 +0200 Subject: [PATCH] Prevent ambiguity --- ext/VoronoiFVMExtendableFEMBaseExt.jl | 63 +++++++++++++-------------- 1 file changed, 30 insertions(+), 33 deletions(-) diff --git a/ext/VoronoiFVMExtendableFEMBaseExt.jl b/ext/VoronoiFVMExtendableFEMBaseExt.jl index 4ac69f391..f3df4e6fc 100644 --- a/ext/VoronoiFVMExtendableFEMBaseExt.jl +++ b/ext/VoronoiFVMExtendableFEMBaseExt.jl @@ -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 @@ -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) @@ -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 @@ -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 @@ -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) @@ -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, @@ -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 @@ -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