Skip to content

Commit

Permalink
Merge pull request #152 from GalerkinToolkit/integration
Browse files Browse the repository at this point in the history 
save integration update and fix bugs
  • Loading branch information
@fverdugo
fverdugo authored Dec 3, 2024
2 parents 6b072bc + 9788dbf commit 405bac0
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 62 deletions.
9 changes: 8 additions & 1 deletion src/domain.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
101 changes: 47 additions & 54 deletions src/integration.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down
41 changes: 34 additions & 7 deletions src/symbolics.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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...)

Expand Down Expand Up @@ -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...)

Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 405bac0

Please sign in to comment.