From 9788dbf9ebbbd3295a0c7ab8a6c1500f845943c8 Mon Sep 17 00:00:00 2001 From: Xiaowei Date: Tue, 3 Dec 2024 14:27:47 +0100 Subject: [PATCH] save integration --- src/domain.jl | 9 +++- src/integration.jl | 101 +++++++++++++++++++++------------------------ src/symbolics.jl | 41 ++++++++++++++---- 3 files changed, 89 insertions(+), 62 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index 2b459ca..79f3652 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -948,7 +948,14 @@ function physical_map(mesh::AbstractMesh,d) end end -function inverse_physical_map(mesh::AbstractMesh,d) +function physical_map(mesh::AbstractDomain{<:PMesh},d) + quantity() do index + dval = Val(val_parameter(d)) + map(x -> GT.physical_map_term(dval, index),partition(mesh)) + end +end + +function inverse_physical_map(mesh::PMesh,d) x0 = zero(SVector{val_parameter(d),Float64}) x = constant_quantity(x0) phi = physical_map(mesh,d) diff --git a/src/integration.jl b/src/integration.jl index 280639b..5a60a5f 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -185,9 +185,7 @@ end function coordinates(measure::Measure{<:PMesh}) q = map(GT.coordinates,partition(measure)) term = map(GT.term,q) - prototype = map(GT.prototype,q) |> PartitionedArrays.getany - domain = measure.domain - GT.quantity(term,prototype,domain) + GT.quantity(term) end function coordinates(measure::Measure,domain::ReferenceDomain) @@ -198,24 +196,14 @@ function coordinates(measure::Measure,domain::ReferenceDomain) refid_to_quad = reference_quadratures(measure) refid_to_coords = map(GT.coordinates,refid_to_quad) prototype = first(GT.coordinates(first(refid_to_quad))) - GT.quantity(prototype,domain) do index - domface = index.face - point = index.point - domface_to_face_sym = get!(index.dict,domface_to_face,gensym("domface_to_face")) - face_to_refid_sym = get!(index.dict,face_to_refid,gensym("face_to_refid")) - refid_to_coords_sym = get!(index.dict,refid_to_coords,gensym("refid_to_coords")) - @term begin - face = $domface_to_face_sym[$domface] - coords = reference_value($refid_to_coords_sym,$face_to_refid_sym,face) - coords[$point] - end - end + + point_quantity(refid_to_coords,domain;reference=true) end function coordinates(measure::Measure,domain::PhysicalDomain) Ω = domain Ωref = GT.reference_domain(Ω) - ϕ = GT.domain_map(Ωref,Ω) + ϕ = GT.physical_map(mesh(Ω),face_dim(Ω)) x = GT.coordinates(measure,Ωref) ϕ(x) end @@ -241,17 +229,18 @@ function weights(measure::Measure,domain::ReferenceDomain) refid_to_quad = reference_quadratures(measure) refid_to_ws = map(GT.weights,refid_to_quad) prototype = first(GT.weights(first(refid_to_quad))) - GT.quantity(prototype,domain) do index - domface = index.face - point = index.point - domface_to_face_sym = get!(index.dict,domface_to_face,gensym("domface_to_face")) - face_to_refid_sym = get!(index.dict,face_to_refid,gensym("face_to_refid")) - refid_to_ws_sym = get!(index.dict,refid_to_ws,gensym("refid_to_ws")) - @term begin + GT.quantity() do index + domface = face_index(index,d) + point = point_index(index) + domface_to_face_sym = get_symbol!(index,domface_to_face,gensym("domface_to_face")) + face_to_refid_sym = get_symbol!(index,face_to_refid,gensym("face_to_refid")) + refid_to_ws_sym = get_symbol!(index,refid_to_ws,gensym("refid_to_ws")) + expr = @term begin face = $domface_to_face_sym[$domface] - ws = reference_value($refid_to_ws_sym,$face_to_refid_sym,face) + ws = $refid_to_ws_sym[$face_to_refid_sym[face]] ws[$point] end + expr_term(d,expr,prototype,index) end end @@ -262,7 +251,7 @@ end function weights(measure::Measure,domain::PhysicalDomain) Ω = domain Ωref = GT.reference_domain(Ω) - ϕ = GT.domain_map(Ωref,Ω) + ϕ = GT.physical_map(mesh(Ω),face_dim(Ω)) w = GT.weights(measure,Ωref) x = GT.coordinates(measure,Ωref) J = ForwardDiff.jacobian(ϕ,x) @@ -280,15 +269,16 @@ function num_points(measure::Measure) refid_to_quad = reference_quadratures(measure) refid_to_ws = map(GT.weights,refid_to_quad) index -> begin - domface = index.face - domface_to_face_sym = get!(index.dict,domface_to_face,gensym("domface_to_face")) - face_to_refid_sym = get!(index.dict,face_to_refid,gensym("face_to_refid")) - refid_to_ws_sym = get!(index.dict,refid_to_ws,gensym("refid_to_ws")) - @term begin + domface = face_index(index, d) + domface_to_face_sym = get_symbol!(index,domface_to_face,gensym("domface_to_face")) + face_to_refid_sym = get_symbol!(index,face_to_refid,gensym("face_to_refid")) + refid_to_ws_sym = get_symbol!(index,refid_to_ws,gensym("refid_to_ws")) + expr = @term begin face = $domface_to_face_sym[$domface] - ws = reference_value($refid_to_ws_sym,$face_to_refid_sym,face) + ws = $refid_to_ws_sym[$face_to_refid_sym[face]] length(ws) end + expr_term(d,expr,0,index) end end @@ -319,9 +309,7 @@ function integrate(f,measure::Measure{<:PMesh}) fx = f(x) q = map(GT.integrate_impl,partition(fx),partition(measure)) term = map(GT.term,q) - prototype = map(GT.prototype,q) |> PartitionedArrays.getany - domain = measure |> GT.domain - contrib = GT.quantity(term,prototype,domain) + contrib = GT.quantity(term) integral((measure=>contrib,)) end @@ -366,20 +354,22 @@ end function sum_contribution_impl(qty,measure,facemask) # TODO some code duplication with face_contribution_impl # Yes, but this allocates less memory - term_qty = GT.term(qty) - term_npoints = GT.num_points(measure) - face = :face - point = :point - index = GT.index(;face,point) - expr_qty = term_qty(index) |> simplify - expr_npoints = term(npoints,index) |> expression |> simplify + term_npoints = GT.num_points(measure) # TODO: term? + dom = domain(measure) + d = GT.num_dims(dom) + index = GT.generate_index(dom) + t = term(qty, index) + face = face_index(index,d) + point = point_index(index) + expr_qty = t |> expression |> simplify + expr_npoints = term_npoints(index) |> expression |> simplify # TODO merge statements s_qty = GT.topological_sort(expr_qty,(face,point)) s_npoints = GT.topological_sort(expr_npoints,(face,)) expr = quote function loop(z,facemask,storage) s = zero(z) - $(unpack_storage(index.dict,:storage)) + $(unpack_index_storage(index,:storage)) $(s_qty[1]) $(s_npoints[1]) nfaces = length(facemask) @@ -397,8 +387,8 @@ function sum_contribution_impl(qty,measure,facemask) end end loop = eval(expr) - storage = GT.storage(index) - z = zero(GT.prototype(qty)) + storage = GT.index_storage(index) + z = zero(GT.prototype(t)) s = invokelatest(loop,z,facemask,storage) s end @@ -421,19 +411,22 @@ function face_contribution(measure,qty) end function face_contribution_impl(qty,measure,facemask) - term_qty = GT.term(qty) term_npoints = GT.num_points(measure) - face = :face - point = :point - index = GT.index(;face,point) - expr_qty = term_qty(index) |> simplify - expr_npoints = term_npoints(index) |> simplify + dom = domain(measure) + d = GT.num_dims(dom) + index = GT.generate_index(dom) + t = term(qty, index) + face = face_index(index,d) + point = point_index(index) + + expr_qty = t |> expression |> simplify + expr_npoints = term_npoints(index) |> expression |> simplify # TODO merge statements s_qty = GT.topological_sort(expr_qty,(face,point)) s_npoints = GT.topological_sort(expr_npoints,(face,)) expr = quote function loop!(facevals,facemask,storage) - $(unpack_storage(index.dict,:storage)) + $(unpack_index_storage(index,:storage)) $(s_qty[1]) $(s_npoints[1]) nfaces = length(facemask) @@ -454,8 +447,8 @@ function face_contribution_impl(qty,measure,facemask) end end loop! = eval(expr) - storage = GT.storage(index) - z = zero(GT.prototype(qty)) + storage = GT.index_storage(index) + z = zero(GT.prototype(t)) nfaces = length(facemask) facevals = fill(z,nfaces) invokelatest(loop!,facevals,facemask,storage) @@ -465,7 +458,7 @@ end function face_diameter(Ω) dΩ = GT.measure(Ω,1) d = GT.num_dims(Ω) - u = GT.analytical_field(x->1,Ω) + u = GT.analytical_field(x->1) int = ∫(u,dΩ) r = face_contribution(int,Ω) r .= r .^ (1/d) diff --git a/src/symbolics.jl b/src/symbolics.jl index dcd7667..da9de1b 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -162,7 +162,7 @@ end free_dims(a::ConstantTerm) = Int[] -prototype(t::ConstantTerm) = a.value +prototype(a::ConstantTerm) = a.value reference_face_term(args...) = ReferenceFaceTerm(args...) @@ -200,12 +200,35 @@ reference_face_term(args...) = ReferenceFaceTerm(args...) # zero(eltype(a.sface_to_value)) #end -#unary_call_term(args...) = UnaryCallTerm(args...) -# -#struct UnaryCallTerm{A,B} <: AbstractTerm -# callee::A -# arg::B -#end +unary_call_term(args...) = UnaryCallTerm(args...) + +struct UnaryCallTerm{A,B} <: AbstractTerm + callee::A + arg::B +end + +AbstractTrees.children(a::UnaryCallTerm) = (a.callee,a.arg) + +index(a::UnaryCallTerm) = index(a.arg) +function prototype(a::UnaryCallTerm) + return_prototype(a.callee,prototype(a.arg)) +end + +free_dims(a::UnaryCallTerm) = free_dims(a.arg) + + +function expression(c::UnaryCallTerm) + f = c.callee + a = c.arg + arg = expression(a) + # ndofs = expression(a.ndofs) + # dof = a.dof + expr = @term begin + $f($arg) + # fun = $dof -> $c*$s + # sum(fun,1:$ndofs) + end +end binary_call_term(args...) = BinaryCallTerm(args...) @@ -652,6 +675,10 @@ end free_dims(a::PhysicalMapTerm) = [val_parameter(a.dim)] + +prototype(a::PhysicalMapTerm) = x-> zero(SVector{val_parameter(a.dim),Float64}) + + function Base.:(==)(a::PhysicalMapTerm,b::PhysicalMapTerm) a.dim == b.dim end