@@ -39,6 +39,10 @@ struct Polynomial{N, T}
3939 order2coeff:: Dict{NTuple{N, Int}, T}
4040end
4141
42+ function order2coeff (p:: Polynomial{N, T} ) where {N, T}
43+ return p. order2coeff
44+ end
45+
4246# empty constructor
4347Polynomial {N, T} () where {N, T} = Polynomial {N, T} (Dict {NTuple{N, Int}, T} ())
4448
@@ -509,6 +513,7 @@ function solve_laplace(Q::Polynomial{N, T}) where {N, T}
509513 @assert is_homogeneous (Q) " source term `Q` must be a homogeneous polynomial"
510514 n = degree (Q)
511515 γ = (k, p) -> 2 * (k + 1 ) * (2 k + 2 p + N) # γₖᵖ
516+ # Note: convert to big for the intermediate computations, then back to T at the end
512517 cₖ = big (1 ) // γ (0 , n) # c₀
513518 P = cₖ * multiply_by_r (deepcopy (Q), 2 )
514519 ΔᵏQ = deepcopy (Q)
@@ -519,7 +524,7 @@ function solve_laplace(Q::Polynomial{N, T}) where {N, T}
519524 ΔP = cₖ * (multiply_by_r (ΔᵏQ, 2 k + 2 ))
520525 P = P + ΔP
521526 end
522- return P
527+ return convert_coefs (P, T)
523528end
524529
525530"""
@@ -558,7 +563,7 @@ function solve_anisotropic_laplace(A::AbstractMatrix{T}, Q::Polynomial{N, T}) wh
558563 ΔP = cₖ * (multiply_by_anisotropic_r (A, ΔᵏQ, 2 k + 2 ))
559564 P = P + ΔP
560565 end
561- return P
566+ return convert_coefs (P, T)
562567end
563568
564569"""
@@ -628,7 +633,7 @@ function solve_anisotropic_advect(β::AbstractVector, Q::Polynomial{N, T}) where
628633 betagradellq = sum (β[i] * gradient (betagradellq)[i] for i in 1 : N)
629634 P = P + cₗ * multiply_by_anisotropic_β_r (β, betagradellq, l + 1 )
630635 end
631- return (1 / β2) * P
636+ return (1 / β2) * convert_coefs (P, T)
632637end
633638
634639"""
0 commit comments