From 18151d0f0de576138da7107bacb8619cf2899e10 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Mon, 1 Jul 2024 18:15:46 +0200 Subject: [PATCH 01/20] Starting to add code generation [skip ci] --- src/GalerkinToolkit.jl | 1 + src/domain.jl | 193 +++++++++++++++++++++++++--------------- src/symbolics.jl | 62 +++++++++++++ test/domain_tests.jl | 4 +- test/runtests.jl | 1 + test/symbolics_tests.jl | 3 + 6 files changed, 189 insertions(+), 75 deletions(-) create mode 100644 src/symbolics.jl create mode 100644 test/symbolics_tests.jl diff --git a/src/GalerkinToolkit.jl b/src/GalerkinToolkit.jl index 838c6d3b..466630c2 100644 --- a/src/GalerkinToolkit.jl +++ b/src/GalerkinToolkit.jl @@ -14,6 +14,7 @@ using SparseArrays using Metis include("geometry.jl") +include("symbolics.jl") include("domain.jl") include("integration.jl") include("interpolation.jl") diff --git a/src/domain.jl b/src/domain.jl index d72a242d..9880eed0 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -444,19 +444,39 @@ function call(f,args...) f(args...) end -function call(g,args::AbstractQuantity...) +function call(g::AbstractQuantity,args::AbstractQuantity...) fs = map(gk.term,args) domain = args |> first |> gk.domain #msg = "All quantities need to be defined on the same domain" #@assert all(dom->dom==domain,map(gk.domain,args)) msg # TODO check everything except reference/physical domain? # Maybe check reference/physical domain only when evaluating functions? - prototype = gk.return_prototype(g,map(gk.prototype,args)...) + prototype = gk.return_prototype(gk.prototype(g),(map(gk.prototype,args)...)) + g_term = gk.term(g) gk.quantity(prototype,domain) do index - g(map(f->f(index),fs)...) + f_exprs = map(f->f(index),fs) + g_expr = g_term(index) + @term g_expr(f_exprs...) end end +function call(g,args::AbstractQuantity...) + domain = args |> first |> gk.domain + g_qty = gk.constant_quantity(g,domain) + call(g_qty,args...) + #fs = map(gk.term,args) + #domain = args |> first |> gk.domain + ##msg = "All quantities need to be defined on the same domain" + ##@assert all(dom->dom==domain,map(gk.domain,args)) msg + ## TODO check everything except reference/physical domain? + ## Maybe check reference/physical domain only when evaluating functions? + #prototype = gk.return_prototype(g,map(gk.prototype,args)...) + #gk.quantity(prototype,domain) do index + # f_exprs = map(f->f(index),fs) + # @term g(f_exprs...) + #end +end + function call(g,args::AbstractQuantity{<:PMesh}...) pargs = map(partition,args) q = map(pargs...) do myargs... @@ -473,7 +493,7 @@ function (f::AbstractQuantity)(x::AbstractQuantity) codomain = gk.domain(f) flag = physical_domain(domain) == physical_domain(codomain) if flag - call(call,f,x) + call(f,x) else align_and_call(f,x) end @@ -488,12 +508,12 @@ end function align_and_call(f,x,glue::InteriorGlue) g = gk.align_field(f,gk.domain(glue)) - call(call,g,x) + call(g,x) end function align_and_call(f,x,glue::BoundaryGlue) g = gk.align_field(f,gk.domain(glue)) - call(call,g,x) + call(g,x) end function align_and_call(f,x,glue::CoboundaryGlue) @@ -510,7 +530,8 @@ function face_constant_field(data,dom::AbstractDomain) prototype = x->zero(eltype(data)) gk.quantity(prototype,dom) do index face = index.face - x->data[face] + face_constant_call(data,face) = x->data[face] + @term face_constant_call(data,face) end end @@ -540,6 +561,34 @@ function domain_map(glue::InteriorGlue,::PhysicalDomain,::PhysicalDomain) gk.quantity(term,prototype,domain) end +function linear_combination(funs,coefs) + q -> begin + sum(1:length(coefs)) do i + x = coefs[i] + fun = funs[i] + coeff = fun(q) + coeff*x + end + end +end + +function linear_combination(funs,coefs,ids) + q -> begin + sum(1:length(ids)) do i + x = coefs[ids[i]] + fun = funs[i] + coeff = fun(q) + coeff*x + end + end +end + +function reference_value(rid_to_value,face_to_rid,face) + rid = face_to_rid[face] + value = rid_to_value[rid] + value +end + function domain_map(glue::InteriorGlue,::ReferenceDomain,::PhysicalDomain) domain = glue.domain mesh = domain |> gk.mesh @@ -555,18 +604,11 @@ function domain_map(glue::InteriorGlue,::ReferenceDomain,::PhysicalDomain) prototype = y->x gk.quantity(prototype,domain) do index sface = index.face - face = sface_to_face[sface] - refid = face_to_refid[face] - funs = refid_to_funs[refid] - nodes = face_to_nodes[face] - coords = node_to_coords[nodes] - q -> begin - sum(1:length(coords)) do i - x = coords[i] - fun = funs[i] - coeff = fun(q) - coeff*x - end + @term begin + face = sface_to_face[sface] + funs = reference_value(refid_to_funs,face_to_refid,face) + nodes = face_to_nodes[face] + linear_combination(funs,node_to_coords,nodes) end end end @@ -621,27 +663,22 @@ function domain_map(glue::CoboundaryGlue,::ReferenceDomain,::ReferenceDomain) drefid_to_funs = map(gk.shape_functions,drefid_to_refdface) gk.quantity(prototype,domain) do index sface = index.face - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - faces_around = sface_to_faces_around[sface] - n_faces_around = length(lfaces) - map(faces_around) do face_around - tface = tfaces[face_around] - lface = lfaces[face_around] - Dface = tface_to_Dface[tface] - dface = sface_to_dface[sface] - perm = Dface_to_lface_to_perm[Dface][lface] - Drefid = Dface_to_Drefid[Dface] - drefid = dface_to_drefid[dface] - coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] - funs = drefid_to_funs[drefid] - q -> begin - sum(1:length(coords)) do i - x = coords[i] - fun = funs[i] - coeff = fun(q) - coeff*x - end + @term begin + tfaces = sface_to_tfaces[sface] + lfaces = sface_to_lfaces[sface] + faces_around = sface_to_faces_around[sface] + n_faces_around = length(lfaces) + map(faces_around) do face_around + tface = tfaces[face_around] + lface = lfaces[face_around] + Dface = tface_to_Dface[tface] + dface = sface_to_dface[sface] + perm = Dface_to_lface_to_perm[Dface][lface] + Drefid = Dface_to_Drefid[Dface] + drefid = dface_to_drefid[dface] + coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] + funs = drefid_to_funs[drefid] + linear_combination(funs,coords) end end end @@ -701,22 +738,17 @@ function domain_map(glue::BoundaryGlue,::ReferenceDomain,::ReferenceDomain) drefid_to_funs = map(gk.shape_functions,drefid_to_refdface) gk.quantity(prototype,domain) do index sface = index.face - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] - Dface = tface_to_Dface[tface] - dface = sface_to_dface[sface] - perm = Dface_to_lface_to_perm[Dface][lface] - Drefid = Dface_to_Drefid[Dface] - drefid = dface_to_drefid[dface] - coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] - funs = drefid_to_funs[drefid] - q -> begin - sum(1:length(coords)) do i - x = coords[i] - fun = funs[i] - coeff = fun(q) - coeff*x - end + @term begin + tface = sface_to_tfaces[sface][1] + lface = sface_to_lfaces[sface][1] + Dface = tface_to_Dface[tface] + dface = sface_to_dface[sface] + perm = Dface_to_lface_to_perm[Dface][lface] + Drefid = Dface_to_Drefid[Dface] + drefid = dface_to_drefid[dface] + coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] + funs = drefid_to_funs[drefid] + linear_combination(funs,coords) end end end @@ -748,7 +780,7 @@ function align_field(a::AbstractQuantity,glue::InteriorGlue) sface_to_tfaces, = gk.target_face(glue) gk.quantity(prototype,domain) do index sface = index.face - tface = sface_to_tfaces[sface][1] + @term tface = sface_to_tfaces[sface][1] index2 = replace_face(index,tface) ai = term_a(index2) ai @@ -761,6 +793,7 @@ function align_field(a::AbstractQuantity,glue::CoboundaryGlue) domain = glue |> gk.domain term_a = gk.term(a) sface_to_tfaces, sface_to_lfaces, = glue |> gk.target_face + # TODO use @term somewhere gk.quantity(prototype,domain) do index sface = index.face tfaces = sface_to_tfaces[sface] @@ -786,8 +819,8 @@ function align_field(a::AbstractQuantity,glue::BoundaryGlue) face_around = gk.face_around(glue.domain) gk.quantity(prototype,domain) do index sface = index.face - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] + @term tface = sface_to_tfaces[sface][1] + @term lface = sface_to_lfaces[sface][1] index2 = replace_face(index,tface) index3 = replace_face_around(index2,face_around) ai = term_a(index3) @@ -853,11 +886,11 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::InteriorGlue) sface_to_tfaces, = gk.target_face(glue) gk.quantity(prototype,domain) do index sface = index.face - tface = sface_to_tfaces[sface][1] + @term tface = sface_to_tfaces[sface][1] index2 = replace_face(index,tface) ai = term_a(index2) phii = term_phi(index) - x-> ai(phii(x)) + @term ai∘phii end end @@ -871,6 +904,7 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::CoboundaryGlue) term_a = gk.term(a) term_phi = gk.term(phi) sface_to_tfaces, sface_to_lfaces, = glue |> gk.target_face + # TODO use @term somewhere gk.quantity(prototype,domain) do index sface = index.face tfaces = sface_to_tfaces[sface] @@ -901,13 +935,13 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::BoundaryGlue) face_around = gk.face_around(glue.domain) gk.quantity(prototype,domain) do index sface = index.face - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] + @term tface = sface_to_tfaces[sface][1] + @term lface = sface_to_lfaces[sface][1] index2 = replace_face(index,tface) index3 = replace_face_around(index2,face_around) ai = term_a(index3) phii = term_phi(index) - x -> ai(phii(x)) + @term ai∘phii end end @@ -956,11 +990,13 @@ function reference_coordinates(plt::Plot) face_to_refid = gk.face_reference_id(mesh,d) prototype = first(first(refid_to_snode_to_coords)) gk.quantity(prototype,domain) do index - domface = index.face[1] + domface = index.face point = index.point - face = domface_to_face[domface] - refid = face_to_refid[face] - refid_to_snode_to_coords[refid][point] + @term begin + face = domface_to_face[domface] + refid = face_to_refid[face] + refid_to_snode_to_coords[refid][point] + end end end @@ -1005,14 +1041,23 @@ function plot_impl!(plt,term,label,::Type{T}) where T nnodes = gk.num_nodes(vmesh) data = zeros(T,nnodes) face_to_nodes = vglue.face_fine_nodes - for face in 1:length(face_to_nodes) - nodes = face_to_nodes[face] - for point in 1:length(nodes) - index = gk.index(;face,point) - v = term(index) - data[nodes[point]] = v + face = gensym("face") + point = gensym("point") + index = gk.index(;face,point) + v = term(index) |> gk.topological_sort + expr = quote + function filldata!(data,face_to_nodes) + for $face in 1:length(face_to_nodes) + nodes = face_to_nodes[$face] + for $point in 1:length(nodes) + data[nodes[$point]] = $v + end + end end end + display(expr) + filldata! = eval(expr) + invokelatest(filldata!,data,face_to_nodes) plt.node_data[label] = data plt end diff --git a/src/symbolics.jl b/src/symbolics.jl new file mode 100644 index 00000000..803c944c --- /dev/null +++ b/src/symbolics.jl @@ -0,0 +1,62 @@ + +macro term(expr) + transform(a) = a # :($(Expr(:quote,a))) + function transform(expr::Expr) + if expr.head === :call + transform_call(expr) + elseif expr.head === :ref + transform_ref(expr) + else + transform_default(expr) + end + end + function transform_call(expr::Expr) + args = map(transform,expr.args) + quote + Expr(:call,$(args...)) + end + end + function transform_ref(expr::Expr) + args = map(transform,expr.args) + quote + Expr(:ref,$(args...)) + end + end + function transform_default(expr::Expr) + head = expr.head + args = map(transform,expr.args) + Expr(head,args...) + end + transform(expr) |> esc +end + +function topological_sort(expr) + temporary = gensym() + expr_L = Expr(:block) + marks = Dict{UInt,Symbol}() + function visit(expr_n) + id_n = objectid(expr_n) + if haskey(marks,id_n) + if marks[id_n] !== temporary + return marks[id_n] + else + error("Graph has cycles! This is not possible.") + end + end + marks[id_n] = temporary + var = gensym() + if isa(expr_n,Expr) + args = expr_n.args + args_var = map(visit,args) + expr_n_new = Expr(expr_n.head,args_var...) + assignment = :($var = $expr_n_new) + else + assignment = :($var = $expr_n) + end + marks[id_n] = var + push!(expr_L.args,assignment) + var + end + visit(expr) + expr_L +end diff --git a/test/domain_tests.jl b/test/domain_tests.jl index 13025afe..56d58d9b 100644 --- a/test/domain_tests.jl +++ b/test/domain_tests.jl @@ -27,9 +27,11 @@ u2 = uref∘ϕinv gk.vtk_plot(joinpath(outdir,"omega"),Ω;refinement=4) do plt gk.plot!(plt,u;label="u") - gk.plot!(plt,u2;label="u2") + #gk.plot!(plt,u2;label="u2") end +xxxx + gk.vtk_plot(joinpath(outdir,"omega_ref"),Ωref;refinement=4) do plt gk.plot!(plt,uref;label="u") diff --git a/test/runtests.jl b/test/runtests.jl index bac7c2a2..0df26b90 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using Test @testset "GalerkinToolkit" begin @testset "Geometry" begin include("geometry_tests.jl") end + @testset "Symbolics" begin include("symbolics_tests.jl") end @testset "Domain" begin include("domain_tests.jl") end @testset "Integration" begin include("integration_tests.jl") end @testset "Interpolation" begin include("interpolation_tests.jl") end diff --git a/test/symbolics_tests.jl b/test/symbolics_tests.jl new file mode 100644 index 00000000..2d8c7675 --- /dev/null +++ b/test/symbolics_tests.jl @@ -0,0 +1,3 @@ +module SymbolicsTests + +end # module From 4410d139a5937a09166515d59b15f97dd41e41f6 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 3 Jul 2024 14:11:17 +0200 Subject: [PATCH 02/20] Implemented some vanilla loop hoisting --- src/domain.jl | 7 +++++-- src/symbolics.jl | 25 +++++++++++++++++++------ 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index 9880eed0..80cb2999 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -1044,13 +1044,16 @@ function plot_impl!(plt,term,label,::Type{T}) where T face = gensym("face") point = gensym("point") index = gk.index(;face,point) - v = term(index) |> gk.topological_sort + deps = (face,point) + v = gk.topological_sort(term(index),deps) expr = quote function filldata!(data,face_to_nodes) + $(v[1]) for $face in 1:length(face_to_nodes) nodes = face_to_nodes[$face] + $(v[2]) for $point in 1:length(nodes) - data[nodes[$point]] = $v + data[nodes[$point]] = $(v[3]) end end end diff --git a/src/symbolics.jl b/src/symbolics.jl index 803c944c..83d84f14 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -30,15 +30,24 @@ macro term(expr) transform(expr) |> esc end -function topological_sort(expr) +function topological_sort(expr,deps) temporary = gensym() - expr_L = Expr(:block) + expr_L = [Expr(:block) for _ in 0:length(deps)] marks = Dict{UInt,Symbol}() + marks_deps = Dict{UInt,Int}() + function visit(expr_n::Symbol) + id_n = objectid(expr_n) + marks[id_n] = expr_n + i = findfirst(e->expr_n===e,deps) + j = i === nothing ? 0 : Int(i) + marks_deps[id_n] = j + expr_n , j + end function visit(expr_n) id_n = objectid(expr_n) if haskey(marks,id_n) if marks[id_n] !== temporary - return marks[id_n] + return marks[id_n], marks_deps[id_n] else error("Graph has cycles! This is not possible.") end @@ -47,15 +56,19 @@ function topological_sort(expr) var = gensym() if isa(expr_n,Expr) args = expr_n.args - args_var = map(visit,args) + r = map(visit,args) + args_var = map(first,r) + j = maximum(map(last,r)) expr_n_new = Expr(expr_n.head,args_var...) assignment = :($var = $expr_n_new) else + j = 0 assignment = :($var = $expr_n) end marks[id_n] = var - push!(expr_L.args,assignment) - var + marks_deps[id_n] = j + push!(expr_L[j+1].args,assignment) + var, j end visit(expr) expr_L From ddf0d642037cbc46c0c5f8e7a7d0f65b2b9b5098 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 3 Jul 2024 17:53:05 +0200 Subject: [PATCH 03/20] Starting to add simplification rules --- Project.toml | 5 +- src/GalerkinToolkit.jl | 1 + src/domain.jl | 126 ++++++++++++++++++++++++++++++----------- src/symbolics.jl | 18 ++++++ test/domain_tests.jl | 4 +- 5 files changed, 119 insertions(+), 35 deletions(-) diff --git a/Project.toml b/Project.toml index bba1d3a0..0ec8a5bb 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +Metatheory = "e9d8d322-4543-424a-9be4-0cc815abe26c" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -18,16 +19,16 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] -FastGaussQuadrature = "1" Combinatorics = "1" +FastGaussQuadrature = "1" ForwardDiff = "0.10" Gmsh = "0.3" IterativeSolvers = "0.9" MPI = "0.20" Metis = "1" PartitionedArrays = "0.4" -WriteVTK = "1" StaticArrays = "1" +WriteVTK = "1" julia = "1" [extras] diff --git a/src/GalerkinToolkit.jl b/src/GalerkinToolkit.jl index 466630c2..23495e0f 100644 --- a/src/GalerkinToolkit.jl +++ b/src/GalerkinToolkit.jl @@ -12,6 +12,7 @@ using PartitionedArrays using Combinatorics using SparseArrays using Metis +using Metatheory include("geometry.jl") include("symbolics.jl") diff --git a/src/domain.jl b/src/domain.jl index 80cb2999..15aebd62 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -322,7 +322,8 @@ end function constant_quantity(v,domain::AbstractDomain) gk.quantity(v,domain) do index - v + v_sym = get!(index.dict,v,gensym("constant_quantity_value")) + v_sym end end @@ -341,7 +342,9 @@ function index(; point=nothing, field_per_dim =nothing, dof_per_dim=nothing, - face_around_per_dim=nothing + face_around_per_dim=nothing, + state=gensym("state"), + dict=IdDict{Any,Symbol}() ) Index( face, @@ -350,7 +353,9 @@ function index(; point, field_per_dim, dof_per_dim, - face_around_per_dim + face_around_per_dim, + state, + dict, ) end @@ -362,6 +367,12 @@ struct Index{A,B,D,E,F,G,H} field_per_dim::F dof_per_dim::G face_around_per_dim::H + state::Symbol + dict::IdDict{Any,Symbol} +end + +function storage(index::Index) + (;( key=>val for (val,key) in index.dict )...) end function replace_face(index::Index,face) @@ -372,7 +383,9 @@ function replace_face(index::Index,face) index.point, index.field_per_dim, index.dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -384,7 +397,9 @@ function replace_local_face(index::Index,local_face) index.point, index.field_per_dim, index.dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -396,7 +411,9 @@ function replace_face_around(index::Index,face_around) index.point, index.field_per_dim, index.dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -408,7 +425,9 @@ function replace_point(index::Index,point) point, index.field_per_dim, index.dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -420,7 +439,9 @@ function replace_field_per_dim(index::Index,field_per_dim) index.point, field_per_dim, index.dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -432,7 +453,9 @@ function replace_dof_per_dim(index::Index,dof_per_dim) index.point, index.field_per_dim, dof_per_dim, - index.face_around_per_dim + index.face_around_per_dim, + index.state, + index.dict, ) end @@ -456,7 +479,7 @@ function call(g::AbstractQuantity,args::AbstractQuantity...) gk.quantity(prototype,domain) do index f_exprs = map(f->f(index),fs) g_expr = g_term(index) - @term g_expr(f_exprs...) + :($g_expr($(f_exprs...))) end end @@ -572,15 +595,37 @@ function linear_combination(funs,coefs) end end -function linear_combination(funs,coefs,ids) - q -> begin - sum(1:length(ids)) do i - x = coefs[ids[i]] - fun = funs[i] - coeff = fun(q) - coeff*x - end - end +#function linear_combination(funs,coefs,ids) +# q -> begin +# sum(1:length(ids)) do i +# x = coefs[ids[i]] +# fun = funs[i] +# coeff = fun(q) +# coeff*x +# end +# end +#end + +function face_function(rid_to_fs,face_to_rid,face_to_dofs,dofs_to_value,face) + fs = reference_value(rid_to_fs,face_to_rid,face) + dofs = face_to_dofs[face] + values = view(dofs_to_value,dofs) + linear_combination(fs,values) +end + +function reference_point(rid_to_x,face_to_rid,face,point) + reference_value(rid_to_x,face_to_rid,face)[point] +end + +function reference_tabulators(rid_to_fs,rid_to_xs) + map((fs,xs)->map(x->map(f->f(x),fs),xs),rid_to_fs,rid_to_xs) +end + +function face_function_value(rid_to_tab,face_to_rid,face_to_dofs,dofs_to_value,face,point) + tab = reference_value(rid_to_tab,face_to_rid,face) + dofs = face_to_dofs[face] + values = view(dofs_to_value,dofs) + transpose(tab[point])*values end function reference_value(rid_to_value,face_to_rid,face) @@ -604,11 +649,20 @@ function domain_map(glue::InteriorGlue,::ReferenceDomain,::PhysicalDomain) prototype = y->x gk.quantity(prototype,domain) do index sface = index.face - @term begin - face = sface_to_face[sface] - funs = reference_value(refid_to_funs,face_to_refid,face) - nodes = face_to_nodes[face] - linear_combination(funs,node_to_coords,nodes) + state = index.state + dict = index.dict + sface_to_face_sym = get!(dict,sface_to_face,gensym("sface_to_face")) + refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) + face_to_nodes_sym = get!(dict,face_to_nodes,gensym("face_to_nodes")) + node_to_coords_sym = get!(dict,node_to_coords,gensym("node_to_coords")) + face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) + refid_to_refface_sym = get!(dict,refid_to_refface,gensym("face_to_refid")) + @term begin # TODO + face = sface_to_face_sym[sface] + :face_function(refid_to_funs_sym,face_to_refid_sym,face_to_nodes_sym,node_to_coords_sym,face) + #funs = :reference_value(refid_to_funs_sym,face_to_refid_sym,face) + #nodes = face_to_nodes_sym[face] + #:linear_combination(funs,node_to_coords_sym,nodes) end end end @@ -992,10 +1046,14 @@ function reference_coordinates(plt::Plot) gk.quantity(prototype,domain) do index domface = index.face point = index.point + dict = index.dict + domface_to_face_sym = get!(dict,domface_to_face,gensym("domface_to_face")) + face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) + refid_to_snode_to_coords_sym = get!(dict,refid_to_snode_to_coords,gensym("refid_to_snode_to_coords")) @term begin - face = domface_to_face[domface] - refid = face_to_refid[face] - refid_to_snode_to_coords[refid][point] + face = domface_to_face_sym[domface] + coords = :reference_value(refid_to_snode_to_coords_sym,face_to_refid_sym,face) + coords[point] end end end @@ -1036,18 +1094,22 @@ function plot!(plt::Plot,field;label) plot_impl!(plt,term,label,T) end + function plot_impl!(plt,term,label,::Type{T}) where T vmesh,vglue = plt.visualization_mesh nnodes = gk.num_nodes(vmesh) data = zeros(T,nnodes) face_to_nodes = vglue.face_fine_nodes - face = gensym("face") - point = gensym("point") + face = :face + point = :point index = gk.index(;face,point) + index.dict[face_to_nodes] = :face_to_nodes deps = (face,point) - v = gk.topological_sort(term(index),deps) + expr1 = term(index) |> simplify + v = gk.topological_sort(expr1,deps) expr = quote - function filldata!(data,face_to_nodes) + function filldata!(data,state) + (;$(Base.values(index.dict)...)) = state $(v[1]) for $face in 1:length(face_to_nodes) nodes = face_to_nodes[$face] @@ -1060,7 +1122,7 @@ function plot_impl!(plt,term,label,::Type{T}) where T end display(expr) filldata! = eval(expr) - invokelatest(filldata!,data,face_to_nodes) + invokelatest(filldata!,data,gk.storage(index)) plt.node_data[label] = data plt end diff --git a/src/symbolics.jl b/src/symbolics.jl index 83d84f14..76dadf95 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -73,3 +73,21 @@ function topological_sort(expr,deps) visit(expr) expr_L end + +function simplify(expr) + expr +end + +function simplify(expr::Expr) + r001 = @slots a b c d e f g @rule face_function(a,b,c,d,e)(reference_value(f,b,e)[g]) --> face_function_value(reference_tabulators(a,f),b,c,d,e,g) + expr2 = r001(expr) + if expr2 === nothing + args = map(simplify,expr.args) + Expr(expr.head,args...) + else + expr2 + end +end + + + diff --git a/test/domain_tests.jl b/test/domain_tests.jl index 56d58d9b..0118d38d 100644 --- a/test/domain_tests.jl +++ b/test/domain_tests.jl @@ -18,8 +18,10 @@ gk.label_boundary_faces!(mesh;physical_name="boundary_faces") @test Ω == Ω @test Ω != Ωref -u = gk.analytical_field(x->sum(x),Ω) +u = gk.analytical_field(sum,Ω) + ϕ = gk.domain_map(Ωref,Ω) + ϕinv = gk.inverse_map(ϕ) uref = u∘ϕ From 1dd6af93abebb3b59067c9515a448404e1364485 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 3 Jul 2024 18:15:16 +0200 Subject: [PATCH 04/20] Caching gk.faces(domain) --- src/domain.jl | 46 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index 15aebd62..bfc2ea51 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -41,6 +41,7 @@ function domain(mesh; physical_names=gk.physical_names(mesh,face_dim), is_reference_domain = Val(false), face_around=nothing, + cache=nothing, ) if val_parameter(is_reference_domain) @@ -50,6 +51,7 @@ function domain(mesh; physical_names, Val(val_parameter(face_dim)), face_around, + cache, ) else PhysicalDomain( @@ -58,26 +60,41 @@ function domain(mesh; physical_names, Val(val_parameter(face_dim)), face_around, + cache, ) + end |> setup_domain +end + +function setup_domain(domain) + if domain.cache !== nothing + return domain end + faces = gk.faces(domain) + cache = domain_cache(;faces) + replace_cache(domain,cache) +end +function domain_cache(;kwargs...) + (;kwargs...) end -struct PhysicalDomain{A,B,C,D,E} <: AbstractDomain{A} +struct PhysicalDomain{A,B,C,D,E,F} <: AbstractDomain{A} mesh::A mesh_id::B physical_names::C face_dim::Val{D} face_around::E + cache::F end is_reference_domain(a::PhysicalDomain) = false -struct ReferenceDomain{A,B,C,D,E} <: AbstractDomain{A} +struct ReferenceDomain{A,B,C,D,E,F} <: AbstractDomain{A} mesh::A mesh_id::B physical_names::C face_dim::Val{D} face_around::E + cache::F end is_reference_domain(a::ReferenceDomain) = true @@ -96,7 +113,18 @@ function replace_face_around(domain::AbstractDomain,face_around) mesh_id = gk.mesh_id(domain) physical_names = gk.physical_names(domain) is_reference_domain = gk.is_reference_domain(domain) - gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around) + cache = domain.cache + gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around,cache) +end + +function replace_cache(domain::AbstractDomain,cache) + mesh = gk.mesh(domain) + face_dim = gk.face_dim(domain) + mesh_id = gk.mesh_id(domain) + physical_names = gk.physical_names(domain) + is_reference_domain = gk.is_reference_domain(domain) + face_around = gk.face_around(domain) + gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around,cache) end function PartitionedArrays.partition(domain::AbstractDomain{<:PMesh}) @@ -123,7 +151,8 @@ function reference_domain(domain::PhysicalDomain) physical_names = gk.physical_names(domain) is_reference_domain = Val(true) face_around = gk.face_around(domain) - gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around) + cache = domain.cache + gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around,cache) end function reference_domain(domain::ReferenceDomain) @@ -137,7 +166,8 @@ function physical_domain(domain::ReferenceDomain) physical_names = gk.physical_names(domain) face_around = gk.face_around(domain) is_reference_domain = Val(false) - gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around) + cache = domain.cache + gk.domain(mesh;face_dim,mesh_id,physical_names,is_reference_domain,face_around,cache) end function physical_domain(domain::PhysicalDomain) @@ -145,6 +175,9 @@ function physical_domain(domain::PhysicalDomain) end function faces(domain::AbstractDomain) + if domain.cache !== nothing + return domain.cache.faces + end mesh = domain |> gk.mesh D = gk.face_dim(domain) Dface_to_tag = zeros(Int,gk.num_faces(mesh,D)) @@ -660,9 +693,6 @@ function domain_map(glue::InteriorGlue,::ReferenceDomain,::PhysicalDomain) @term begin # TODO face = sface_to_face_sym[sface] :face_function(refid_to_funs_sym,face_to_refid_sym,face_to_nodes_sym,node_to_coords_sym,face) - #funs = :reference_value(refid_to_funs_sym,face_to_refid_sym,face) - #nodes = face_to_nodes_sym[face] - #:linear_combination(funs,node_to_coords_sym,nodes) end end end From 1a156401e691a4b56b2bb3b17bccae268296a6b6 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 3 Jul 2024 18:25:30 +0200 Subject: [PATCH 05/20] Some cosmetic changes in generated code --- src/domain.jl | 5 +++-- src/symbolics.jl | 7 +++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index bfc2ea51..49b7b8ef 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -1133,13 +1133,14 @@ function plot_impl!(plt,term,label,::Type{T}) where T face = :face point = :point index = gk.index(;face,point) - index.dict[face_to_nodes] = :face_to_nodes + dict = index.dict + dict[face_to_nodes] = :face_to_nodes deps = (face,point) expr1 = term(index) |> simplify v = gk.topological_sort(expr1,deps) expr = quote function filldata!(data,state) - (;$(Base.values(index.dict)...)) = state + $(unpack_state(dict,:state)) $(v[1]) for $face in 1:length(face_to_nodes) nodes = face_to_nodes[$face] diff --git a/src/symbolics.jl b/src/symbolics.jl index 76dadf95..0ef7295a 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -89,5 +89,12 @@ function simplify(expr::Expr) end end +function unpack_state(dict,state) + expr = Expr(:block) + for k in Base.values(dict) |> collect |> sort + push!(expr.args,:($k = $state.$k)) + end + expr +end From 9ed30067ab1a4c82e3018a5d9a8b4339ffcc547d Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 3 Jul 2024 19:35:04 +0200 Subject: [PATCH 06/20] Fix type instability --- src/domain.jl | 30 +++++++++++++++++++++++------- src/symbolics.jl | 4 ++-- test/domain_tests.jl | 1 + 3 files changed, 26 insertions(+), 9 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index 49b7b8ef..ca7f8856 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -646,19 +646,35 @@ function face_function(rid_to_fs,face_to_rid,face_to_dofs,dofs_to_value,face) linear_combination(fs,values) end -function reference_point(rid_to_x,face_to_rid,face,point) - reference_value(rid_to_x,face_to_rid,face)[point] -end +#function reference_point(rid_to_x,face_to_rid,face,point) +# reference_value(rid_to_x,face_to_rid,face)[point] +#end + +#function reference_tabulators(rid_to_fs,rid_to_xs) +# map((fs,xs)->map(x->map(f->f(x),fs),xs),rid_to_fs,rid_to_xs) +#end -function reference_tabulators(rid_to_fs,rid_to_xs) - map((fs,xs)->map(x->map(f->f(x),fs),xs),rid_to_fs,rid_to_xs) +function reference_tabulator(fs,xs) + f = fs[1] + z = zero(eltype(xs)) + T = typeof(f(z)) + nx = length(xs) + nf = length(fs) + A = Matrix{T}(undef,nf,nx) + for j in 1:nx + for i in 1:nf + A[i,j] = fs[i](xs[j]) + end + end + A end function face_function_value(rid_to_tab,face_to_rid,face_to_dofs,dofs_to_value,face,point) tab = reference_value(rid_to_tab,face_to_rid,face) dofs = face_to_dofs[face] values = view(dofs_to_value,dofs) - transpose(tab[point])*values + n = length(values) + sum(i->tab[i,point]*values[i],1:n) end function reference_value(rid_to_value,face_to_rid,face) @@ -1140,7 +1156,7 @@ function plot_impl!(plt,term,label,::Type{T}) where T v = gk.topological_sort(expr1,deps) expr = quote function filldata!(data,state) - $(unpack_state(dict,:state)) + $(unpack_storage(dict,:state)) $(v[1]) for $face in 1:length(face_to_nodes) nodes = face_to_nodes[$face] diff --git a/src/symbolics.jl b/src/symbolics.jl index 0ef7295a..b6ac3ec1 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -79,7 +79,7 @@ function simplify(expr) end function simplify(expr::Expr) - r001 = @slots a b c d e f g @rule face_function(a,b,c,d,e)(reference_value(f,b,e)[g]) --> face_function_value(reference_tabulators(a,f),b,c,d,e,g) + r001 = @slots a b c d e f g @rule face_function(a,b,c,d,e)(reference_value(f,b,e)[g]) --> face_function_value(map(reference_tabulator,a,f),b,c,d,e,g) expr2 = r001(expr) if expr2 === nothing args = map(simplify,expr.args) @@ -89,7 +89,7 @@ function simplify(expr::Expr) end end -function unpack_state(dict,state) +function unpack_storage(dict,state) expr = Expr(:block) for k in Base.values(dict) |> collect |> sort push!(expr.args,:($k = $state.$k)) diff --git a/test/domain_tests.jl b/test/domain_tests.jl index 0118d38d..6342b992 100644 --- a/test/domain_tests.jl +++ b/test/domain_tests.jl @@ -32,6 +32,7 @@ gk.vtk_plot(joinpath(outdir,"omega"),Ω;refinement=4) do plt #gk.plot!(plt,u2;label="u2") end + xxxx From b7bf7f540481e7d2abbe4dd92c5425806ae1b5b3 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Sat, 6 Jul 2024 17:47:02 +0200 Subject: [PATCH 07/20] Adding symbolic manipulations to domain.jl --- src/domain.jl | 432 ++++++++++++++++++++++++------------------- src/symbolics.jl | 87 ++++++++- test/domain_tests.jl | 46 ++--- 3 files changed, 345 insertions(+), 220 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index aefbce42..37df0e78 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -400,7 +400,7 @@ struct Index{A,B,D,E,F,G,H} field_per_dim::F dof_per_dim::G face_around_per_dim::H - state::Symbol + state::Symbol #TODO remove dict::IdDict{Any,Symbol} end @@ -516,27 +516,11 @@ function call(g::AbstractQuantity,args::AbstractQuantity...) end end -function call(g,args::AbstractQuantity...) - domain = args |> first |> GT.domain - g_qty = GT.constant_quantity(g,domain) - call(g_qty,args...) - #fs = map(GT.term,args) - #domain = args |> first |> GT.domain - ##msg = "All quantities need to be defined on the same domain" - ##@assert all(dom->dom==domain,map(GT.domain,args)) msg - ## TODO check everything except reference/physical domain? - ## Maybe check reference/physical domain only when evaluating functions? - #prototype = GT.return_prototype(g,map(GT.prototype,args)...) - #GT.quantity(prototype,domain) do index - # f_exprs = map(f->f(index),fs) - # @term g(f_exprs...) - #end -end - -function call(g,args::AbstractQuantity{<:PMesh}...) +function call(g::AbstractQuantity,args::AbstractQuantity{<:PMesh}...) + pg = partition(g) pargs = map(partition,args) - q = map(pargs...) do myargs... - GT.call(g,myargs...) + q = map(pg,pargs...) do myg,myargs... + GT.call(myg,myargs...) end domain = args |> first |> GT.domain term = map(GT.term,q) @@ -544,6 +528,12 @@ function call(g,args::AbstractQuantity{<:PMesh}...) GT.quantity(term,prototype,domain) end +function call(g,args::AbstractQuantity...) + domain = args |> first |> GT.domain + g_qty = GT.constant_quantity(g,domain) + call(g_qty,args...) +end + function (f::AbstractQuantity)(x::AbstractQuantity) domain = GT.domain(x) codomain = GT.domain(f) @@ -573,21 +563,42 @@ function align_and_call(f,x,glue::BoundaryGlue) end function align_and_call(f,x,glue::CoboundaryGlue) - aligned_call(g,x) = map(gi->gi(x),g) g = GT.align_field(f,GT.domain(glue)) - call(aligned_call,g,x) + call(g,x) +end + +function Base.getindex(q::AbstractQuantity,::typeof(+)) + face_around = 1 + term = GT.term(q) + GT.quantity(GT.prototype(q),GT.domain(q)) do index + index2 = GT.replace_face_around(index,face_around) + term(index2) + end +end + +function Base.getindex(q::AbstractQuantity,::typeof(-)) + face_around = 2 + term = GT.term(q) + GT.quantity(GT.prototype(q),GT.domain(q)) do index + index2 = GT.replace_face_around(index,face_around) + term(index2) + end end function analytical_field(f,dom::AbstractDomain) constant_quantity(f,dom) end +function face_constant_field_impl(date,face::Integer) + x->data[face] +end + function face_constant_field(data,dom::AbstractDomain) prototype = x->zero(eltype(data)) GT.quantity(prototype,dom) do index face = index.face - face_constant_call(data,face) = x->data[face] - @term face_constant_call(data,face) + data_sym = get!(index.dict,data,gensym("face_to_constant_value")) + @term face_constant_field_impl($data_sym,$face) end end @@ -606,14 +617,14 @@ end function domain_map(glue::InteriorGlue,::ReferenceDomain,::ReferenceDomain) domain = glue.domain prototype = identity - term = identity + term = index -> :identity GT.quantity(term,prototype,domain) end function domain_map(glue::InteriorGlue,::PhysicalDomain,::PhysicalDomain) domain = glue.domain prototype = identity - term = identity + term = index -> :identity GT.quantity(term,prototype,domain) end @@ -706,9 +717,14 @@ function domain_map(glue::InteriorGlue,::ReferenceDomain,::PhysicalDomain) node_to_coords_sym = get!(dict,node_to_coords,gensym("node_to_coords")) face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) refid_to_refface_sym = get!(dict,refid_to_refface,gensym("face_to_refid")) - @term begin # TODO - face = sface_to_face_sym[sface] - :face_function(refid_to_funs_sym,face_to_refid_sym,face_to_nodes_sym,node_to_coords_sym,face) + @term begin + face = $sface_to_face_sym[$sface] + face_function( + $refid_to_funs_sym, + $face_to_refid_sym, + $face_to_nodes_sym, + $node_to_coords_sym, + $face) end end end @@ -732,8 +748,8 @@ function domain_map(glue::CoboundaryGlue,::ReferenceDomain,::ReferenceDomain) node_to_coords = GT.node_coordinates(boundary) T = eltype(node_to_coords) x = zero(T) - prototype = [y->x,y->x] - sface_to_tfaces, sface_to_lfaces, sface_to_faces_around = glue |> GT.target_face + prototype = y -> x + sface_to_tfaces, sface_to_lfaces, = glue |> GT.target_face tface_to_Dface = codomain |> GT.faces d = domain |> GT.face_dim topo = mesh |> GT.topology @@ -763,23 +779,32 @@ function domain_map(glue::CoboundaryGlue,::ReferenceDomain,::ReferenceDomain) drefid_to_funs = map(GT.shape_functions,drefid_to_refdface) GT.quantity(prototype,domain) do index sface = index.face + dict = index.dict + sface_to_tfaces_sym = get!(dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_Dface_sym = get!(dict,tface_to_Dface,gensym("tface_to_Dface")) + sface_to_dface_sym = get!(dict,sface_to_dface,gensym("sface_to_dface")) + Dface_to_lface_to_perm_sym = get!(dict,Dface_to_lface_to_perm,gensym("Dface_to_lface_to_perm")) + Dface_to_Drefid_sym = get!(dict,Dface_to_Drefid,gensym("Dface_to_Drefid")) + dface_to_drefid_sym = get!(dict,dface_to_drefid,gensym("dface_to_drefid")) + Drefid_to_lface_to_perm_to_coords_sym = get!(dict,Drefid_to_lface_to_perm_to_coords,gensym("Drefid_to_lface_to_perm_to_coords")) + drefid_to_funs_sym = get!(dict,drefid_to_funs,gensym("drefid_to_funs")) + face_around = index.face_around + @assert face_around !== nothing @term begin - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - faces_around = sface_to_faces_around[sface] - n_faces_around = length(lfaces) - map(faces_around) do face_around - tface = tfaces[face_around] - lface = lfaces[face_around] - Dface = tface_to_Dface[tface] - dface = sface_to_dface[sface] - perm = Dface_to_lface_to_perm[Dface][lface] - Drefid = Dface_to_Drefid[Dface] - drefid = dface_to_drefid[dface] - coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] - funs = drefid_to_funs[drefid] - linear_combination(funs,coords) - end + boundary_face_function( + $sface, + $sface_to_tfaces_sym, + $sface_to_lfaces_sym, + $tface_to_Dface_sym, + $sface_to_dface_sym, + $Dface_to_lface_to_perm_sym, + $Dface_to_Drefid_sym, + $dface_to_drefid_sym, + $Drefid_to_lface_to_perm_to_coords_sym, + $drefid_to_funs_sym, + $face_around, + ) end end end @@ -796,6 +821,33 @@ function domain_map(glue::BoundaryGlue,::PhysicalDomain,::PhysicalDomain) error("Case not yet implemented") end +function boundary_face_function( + sface, + sface_to_tfaces, + sface_to_lfaces, + tface_to_Dface, + sface_to_dface, + Dface_to_lface_to_perm, + Dface_to_Drefid, + dface_to_drefid, + Drefid_to_lface_to_perm_to_coords, + drefid_to_funs, + face_around, + ) + tfaces = sface_to_tfaces[sface] + lfaces = sface_to_lfaces[sface] + tface = tfaces[face_around] + lface = lfaces[face_around] + Dface = tface_to_Dface[tface] + dface = sface_to_dface[sface] + perm = Dface_to_lface_to_perm[Dface][lface] + Drefid = Dface_to_Drefid[Dface] + drefid = dface_to_drefid[dface] + coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] + funs = drefid_to_funs[drefid] + linear_combination(funs,coords) +end + function domain_map(glue::BoundaryGlue,::ReferenceDomain,::ReferenceDomain) domain = glue.domain codomain = glue |> GT.codomain @@ -838,17 +890,31 @@ function domain_map(glue::BoundaryGlue,::ReferenceDomain,::ReferenceDomain) drefid_to_funs = map(GT.shape_functions,drefid_to_refdface) GT.quantity(prototype,domain) do index sface = index.face + dict = index.dict + sface_to_tfaces_sym = get!(dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_Dface_sym = get!(dict,tface_to_Dface,gensym("tface_to_Dface")) + sface_to_dface_sym = get!(dict,sface_to_dface,gensym("sface_to_dface")) + Dface_to_lface_to_perm_sym = get!(dict,Dface_to_lface_to_perm,gensym("Dface_to_lface_to_perm")) + Dface_to_Drefid_sym = get!(dict,Dface_to_Drefid,gensym("Dface_to_Drefid")) + dface_to_drefid_sym = get!(dict,dface_to_drefid,gensym("dface_to_drefid")) + Drefid_to_lface_to_perm_to_coords_sym = get!(dict,Drefid_to_lface_to_perm_to_coords,gensym("Drefid_to_lface_to_perm_to_coords")) + drefid_to_funs_sym = get!(dict,drefid_to_funs,gensym("drefid_to_funs")) + face_around = 1 @term begin - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] - Dface = tface_to_Dface[tface] - dface = sface_to_dface[sface] - perm = Dface_to_lface_to_perm[Dface][lface] - Drefid = Dface_to_Drefid[Dface] - drefid = dface_to_drefid[dface] - coords = Drefid_to_lface_to_perm_to_coords[Drefid][lface][perm] - funs = drefid_to_funs[drefid] - linear_combination(funs,coords) + boundary_face_function( + $sface, + $sface_to_tfaces_sym, + $sface_to_lfaces_sym, + $tface_to_Dface_sym, + $sface_to_dface_sym, + $Dface_to_lface_to_perm_sym, + $Dface_to_Drefid_sym, + $dface_to_drefid_sym, + $Drefid_to_lface_to_perm_to_coords_sym, + $drefid_to_funs_sym, + $face_around, + ) end end end @@ -880,7 +946,8 @@ function align_field(a::AbstractQuantity,glue::InteriorGlue) sface_to_tfaces, = GT.target_face(glue) GT.quantity(prototype,domain) do index sface = index.face - @term tface = sface_to_tfaces[sface][1] + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) ai = term_a(index2) ai @@ -888,26 +955,20 @@ function align_field(a::AbstractQuantity,glue::InteriorGlue) end function align_field(a::AbstractQuantity,glue::CoboundaryGlue) - pa = GT.prototype(a) - prototype = [pa,pa] + prototype = GT.prototype(a) domain = glue |> GT.domain term_a = GT.term(a) sface_to_tfaces, sface_to_lfaces, = glue |> GT.target_face - # TODO use @term somewhere GT.quantity(prototype,domain) do index + face_around = index.face_around + @assert face_around !== nothing sface = index.face - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - n_faces_around = length(tfaces) - # TODO This should be a tuple - map(1:n_faces_around) do face_around - tface = sface_to_tfaces[sface][face_around] - lface = sface_to_lfaces[sface][face_around] - index2 = replace_face(index,tface) - index3 = replace_face_around(index2,face_around) - ai = term_a(index3) - ai - end + dict = index.dict + sface_to_tfaces_sym = get!(dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][$face_around] + index2 = replace_face(index,tface) + ai = term_a(index2) + ai end end @@ -919,8 +980,9 @@ function align_field(a::AbstractQuantity,glue::BoundaryGlue) face_around = GT.face_around(glue.domain) GT.quantity(prototype,domain) do index sface = index.face - @term tface = sface_to_tfaces[sface][1] - @term lface = sface_to_lfaces[sface][1] + dict = index.dict + sface_to_tfaces_sym = get!(dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) index3 = replace_face_around(index2,face_around) ai = term_a(index3) @@ -986,11 +1048,12 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::InteriorGlue) sface_to_tfaces, = GT.target_face(glue) GT.quantity(prototype,domain) do index sface = index.face - @term tface = sface_to_tfaces[sface][1] + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) ai = term_a(index2) phii = term_phi(index) - @term ai∘phii + @term $ai∘$phii end end @@ -999,27 +1062,22 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::CoboundaryGlue) @assert GT.domain(a) == GT.codomain(glue) g = GT.prototype(a) f = GT.prototype(phi) - prototype = map(fi->(x->g(fi(x))),f) + prototype = x-> g(f(x)) domain = phi |> GT.domain term_a = GT.term(a) term_phi = GT.term(phi) sface_to_tfaces, sface_to_lfaces, = glue |> GT.target_face - # TODO use @term somewhere + face_around = GT.face_around(glue.domain) GT.quantity(prototype,domain) do index + face_around = index.face_around + @assert face_around !== nothing sface = index.face - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - n_faces_around = length(tfaces) + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][$face_around] + index2 = replace_face(index,tface) + ai = term_a(index2) phii = term_phi(index) - # TODO This should be a tuple - map(1:n_faces_around) do face_around - tface = sface_to_tfaces[sface][face_around] - lface = sface_to_lfaces[sface][face_around] - index2 = replace_face(index,tface) - index3 = replace_face_around(index2,face_around) - ai = term_a(index3) - x -> ai(phii[face_around](x)) - end + @term $ai∘$phii end end @@ -1035,13 +1093,13 @@ function compose(a::AbstractQuantity,phi::AbstractQuantity,glue::BoundaryGlue) face_around = GT.face_around(glue.domain) GT.quantity(prototype,domain) do index sface = index.face - @term tface = sface_to_tfaces[sface][1] - @term lface = sface_to_lfaces[sface][1] + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) index3 = replace_face_around(index2,face_around) ai = term_a(index3) phii = term_phi(index) - @term ai∘phii + @term $ai∘$phii end end @@ -1097,9 +1155,12 @@ function reference_coordinates(plt::Plot) face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) refid_to_snode_to_coords_sym = get!(dict,refid_to_snode_to_coords,gensym("refid_to_snode_to_coords")) @term begin - face = domface_to_face_sym[domface] - coords = :reference_value(refid_to_snode_to_coords_sym,face_to_refid_sym,face) - coords[point] + face = $domface_to_face_sym[$domface] + coords = reference_value( + $refid_to_snode_to_coords_sym, + $face_to_refid_sym, + $face) + coords[$point] end end end @@ -1167,7 +1228,6 @@ function plot_impl!(plt,term,label,::Type{T}) where T end end end - display(expr) filldata! = eval(expr) invokelatest(filldata!,data,GT.storage(index)) plt.node_data[label] = data @@ -1246,7 +1306,23 @@ function unit_normal(domain::AbstractDomain,codomain::AbstractDomain) unit_normal(domain,codomain,glue) end -# TODO a lot of code duplication +# TODO think about this name +function boundary_unit_normal(φ_fun,ϕ_fun,n) + q -> begin + p = φ_fun(q) + J = ForwardDiff.jacobian(ϕ_fun,p) + Jt = transpose(J) + pinvJt = transpose(inv(Jt*J)*Jt) + v = pinvJt*n + m = sqrt(inner(v,v)) + if m < eps() + return zero(v) + else + return v/m + end + end +end + function unit_normal(domain::ReferenceDomain,codomain::PhysicalDomain,glue::BoundaryGlue) Γref = domain Ω = codomain @@ -1268,27 +1344,22 @@ function unit_normal(domain::ReferenceDomain,codomain::PhysicalDomain,glue::Boun prototype = GT.prototype(φ) GT.quantity(prototype,Γref) do index sface = index.face - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] - face = tface_to_face[tface] - ctype = face_to_ctype[face] - lface_to_n = ctype_to_lface_to_n[ctype] - n = lface_to_n[lface] + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(index.dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_face_sym = get!(index.dict,tface_to_face,gensym("tface_to_face")) + face_to_ctype_sym = get!(index.dict,face_to_ctype,gensym("face_to_ctype")) + ctype_to_lface_to_n_sym = get!(index.dict,ctype_to_lface_to_n,gensym("ctype_to_lface_to_n")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) ϕ_fun = ϕ_term(index2) φ_fun = φ_term(index) - q -> begin - p = φ_fun(q) - J = ForwardDiff.jacobian(ϕ_fun,p) - Jt = transpose(J) - pinvJt = transpose(inv(Jt*J)*Jt) - v = pinvJt*n - m = sqrt(inner(v,v)) - if m < eps() - return zero(v) - else - return v/m - end + @term begin + lface = $sface_to_lfaces_sym[$sface][1] + face = $tface_to_face_sym[$tface] + ctype = $face_to_ctype_sym[face] + lface_to_n = $ctype_to_lface_to_n_sym[ctype] + n = lface_to_n[lface] + boundary_unit_normal($φ_fun,$ϕ_fun,n) end end end @@ -1317,27 +1388,22 @@ function unit_normal(domain::PhysicalDomain,codomain::PhysicalDomain,glue::Bound prototype = GT.prototype(ϕ) GT.quantity(prototype,Γref) do index sface = index.face - tface = sface_to_tfaces[sface][1] - lface = sface_to_lfaces[sface][1] - face = tface_to_face[tface] - ctype = face_to_ctype[face] - lface_to_n = ctype_to_lface_to_n[ctype] - n = lface_to_n[lface] + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(index.dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_face_sym = get!(index.dict,tface_to_face,gensym("tface_to_face")) + face_to_ctype_sym = get!(index.dict,face_to_ctype,gensym("face_to_ctype")) + ctype_to_lface_to_n_sym = get!(index.dict,ctype_to_lface_to_n,gensym("ctype_to_lface_to_n")) + tface = @term $sface_to_tfaces_sym[$sface][1] index2 = replace_face(index,tface) ϕinv_fun = ϕinv_term(index2) ϕ_fun = ϕ_term(index2) - x -> begin - q = ϕinv_fun(x) - J = ForwardDiff.jacobian(ϕ_fun,q) - Jt = transpose(J) - pinvJt = transpose(inv(Jt*J)*Jt) - v = pinvJt*n - m = sqrt(inner(v,v)) - if m < eps() - return zero(v) - else - return v/m - end + @term begin + lface = $sface_to_lfaces_sym[$sface][1] + face = $tface_to_face_sym[$tface] + ctype = $face_to_ctype_sym[face] + lface_to_n = $ctype_to_lface_to_n_sym[ctype] + n = lface_to_n[lface] + boundary_unit_normal($ϕinv_fun,$ϕ_fun,n) end end end @@ -1361,44 +1427,32 @@ function unit_normal(domain::ReferenceDomain,codomain::PhysicalDomain,glue::Cobo end ϕ_term = GT.term(ϕ) φ_term = GT.term(φ) - #fun_φ = GT.prototype(φ) - #prototype = [q->(fun_φ(q)[1]),q->(fun_φ(q)[1])] prototype = GT.prototype(φ) GT.quantity(prototype,Γref) do index sface = index.face - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - n_faces_around = length(tfaces) - map(1:n_faces_around) do face_around - tface = tfaces[face_around] - lface = lfaces[face_around] - face = tface_to_face[tface] - ctype = face_to_ctype[face] - lface_to_n = ctype_to_lface_to_n[ctype] + face_around = index.face_around + @assert face_around !== nothing + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(index.dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_face_sym = get!(index.dict,tface_to_face,gensym("tface_to_face")) + face_to_ctype_sym = get!(index.dict,face_to_ctype,gensym("face_to_ctype")) + ctype_to_lface_to_n_sym = get!(index.dict,ctype_to_lface_to_n,gensym("ctype_to_lface_to_n")) + tface = @term $sface_to_tfaces_sym[$sface][$face_around] + index2 = replace_face(index,tface) + index3 = replace_face_around(index2,nothing) + ϕ_fun = ϕ_term(index3) + φ_fun = φ_term(index) + @term begin + lface = $sface_to_lfaces_sym[$sface][$face_around] + face = $tface_to_face_sym[$tface] + ctype = $face_to_ctype_sym[face] + lface_to_n = $ctype_to_lface_to_n_sym[ctype] n = lface_to_n[lface] - index2 = replace_face(index,tface) - ϕ_fun = ϕ_term(index2) - φ_fun = φ_term(index) - # TODO this is not consistent with align_map - # which one is the good one? - q -> begin - p = φ_fun[face_around](q) - J = ForwardDiff.jacobian(ϕ_fun,p) - Jt = transpose(J) - pinvJt = transpose(inv(Jt*J)*Jt) - v = pinvJt*n - m = sqrt(inner(v,v)) - if m < eps() - return zero(v) - else - return v/m - end - end + boundary_unit_normal($φ_fun,$ϕ_fun,n) end end end -# TODO a lot of code duplication function unit_normal(domain::PhysicalDomain,codomain::PhysicalDomain,glue::CoboundaryGlue) Γ = domain Γref = GT.physical_domain(Γ) @@ -1413,40 +1467,34 @@ function unit_normal(domain::PhysicalDomain,codomain::PhysicalDomain,glue::Cobou tface_to_face = GT.faces(Ωref) face_to_ctype = GT.face_reference_id(mesh,D) ctype_to_refface = GT.reference_faces(mesh,D) - ctype_to_lface_to_n = map(ctype_to_refface) do refface + ctype_to_lface_to_n= map(ctype_to_refface) do refface boundary = refface |> GT.geometry |> GT.boundary boundary |> GT.outwards_normals # TODO also rename? end ϕinv_term = GT.term(ϕinv) ϕ_term = GT.term(ϕ) - prototype = [GT.prototype(ϕ),GT.prototype(ϕ)] + prototype = GT.prototype(ϕ) GT.quantity(prototype,Γref) do index sface = index.face - tfaces = sface_to_tfaces[sface] - lfaces = sface_to_lfaces[sface] - map(tfaces,lfaces) do tface,lface - face = tface_to_face[tface] - ctype = face_to_ctype[face] - lface_to_n = ctype_to_lface_to_n[ctype] + face_around = index.face_around + @assert face_around !== nothing + sface_to_tfaces_sym = get!(index.dict,sface_to_tfaces,gensym("sface_to_tfaces")) + sface_to_lfaces_sym = get!(index.dict,sface_to_lfaces,gensym("sface_to_lfaces")) + tface_to_face_sym = get!(index.dict,tface_to_face,gensym("tface_to_face")) + face_to_ctype_sym = get!(index.dict,face_to_ctype,gensym("face_to_ctype")) + ctype_to_lface_to_n_sym = get!(index.dict,ctype_to_lface_to_n,gensym("ctype_to_lface_to_n")) + tface = @term $sface_to_tfaces_sym[$sface][$face_around] + index2 = replace_face(index,tface) + index3 = replace_face_around(index2,nothing) + ϕinv_fun = ϕinv_term(index3) + ϕ_fun = ϕ_term(index3) + @term begin + lface = $sface_to_lfaces_sym[$sface][$face_around] + face = $tface_to_face_sym[$tface] + ctype = $face_to_ctype_sym[face] + lface_to_n = $ctype_to_lface_to_n_sym[ctype] n = lface_to_n[lface] - index2 = replace_face(index,tface) - ϕinv_fun = ϕinv_term(index2) - ϕ_fun = ϕ_term(index2) - # TODO this is not consistent with align_map - # which one is the good one? - x -> begin - q = ϕinv_fun(x) - J = ForwardDiff.jacobian(ϕ_fun,q) - Jt = transpose(J) - pinvJt = transpose(inv(Jt*J)*Jt) - v = pinvJt*n - m = sqrt(inner(v,v)) - if m < eps() - return zero(v) - else - return v/m - end - end + boundary_unit_normal($ϕinv_fun,$ϕ_fun,n) end end end diff --git a/src/symbolics.jl b/src/symbolics.jl index b6ac3ec1..0f4ce0b2 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -1,13 +1,76 @@ +#macro term(expr) +# transform(a) = a # :($(Expr(:quote,a))) +# function transform(expr::Expr) +# if expr.head === :call +# transform_call(expr) +# elseif expr.head === :ref +# transform_ref(expr) +# else +# transform_default(expr) +# end +# end +# function transform_call(expr::Expr) +# args = map(transform,expr.args) +# quote +# Expr(:call,$(args...)) +# end +# end +# function transform_ref(expr::Expr) +# args = map(transform,expr.args) +# quote +# Expr(:ref,$(args...)) +# end +# end +# function transform_default(expr::Expr) +# head = expr.head +# args = map(transform,expr.args) +# Expr(head,args...) +# end +# transform(expr) |> esc +#end + macro term(expr) - transform(a) = a # :($(Expr(:quote,a))) + vars = Set{Symbol}() + function findvars!(a) + nothing + end + function findvars!(expr::Expr) + if expr.head === :(=) + var = expr.args[1] + push!(vars,var) + else + map(findvars!,expr.args) + end + nothing + end + transform(a) = a + function transform(a::Symbol) + if a in vars + a + else + :($(Expr(:quote,a))) + end + end function transform(expr::Expr) if expr.head === :call transform_call(expr) elseif expr.head === :ref transform_ref(expr) - else + elseif expr.head === :$ + transform_interpolation(expr) + elseif expr.head === :(=) transform_default(expr) + elseif expr.head === :(->) + transform_lambda(expr) + elseif expr.head === :do + transform_do(expr) + elseif expr.head === :block + transform_default(expr) + elseif expr.head === :tuple + transform_tuple(expr) + else + error("Expr with head=$(expr.head) not supported in macro @term") end end function transform_call(expr::Expr) @@ -22,11 +85,31 @@ macro term(expr) Expr(:ref,$(args...)) end end + function transform_lambda(expr::Expr) + args = map(transform,expr.args) + quote + Expr(:(->),$(args...)) + end + end + function transform_do(expr::Expr) + expr2 = Expr(:call,expr.args[1].args[1],expr.args[2],expr.args[1].args[2]) + transform(expr2) + end + function transform_interpolation(expr::Expr) + expr.args[1] + end + function transform_tuple(expr::Expr) + args = map(transform,expr.args) + quote + Expr(:tuple,$(args...)) + end + end function transform_default(expr::Expr) head = expr.head args = map(transform,expr.args) Expr(head,args...) end + findvars!(expr) transform(expr) |> esc end diff --git a/test/domain_tests.jl b/test/domain_tests.jl index 4884527c..f4569e93 100644 --- a/test/domain_tests.jl +++ b/test/domain_tests.jl @@ -29,13 +29,9 @@ u2 = uref∘ϕinv GT.vtk_plot(joinpath(outdir,"omega"),Ω;refinement=4) do plt GT.plot!(plt,u;label="u") - #GT.plot!(plt,u2;label="u2") + GT.plot!(plt,u2;label="u2") end - -xxxx - - GT.vtk_plot(joinpath(outdir,"omega_ref"),Ωref;refinement=4) do plt GT.plot!(plt,uref;label="u") end @@ -53,13 +49,13 @@ g = uref∘ϕ n = GT.unit_normal(Γref,Ω) n2 = GT.unit_normal(Γ,Ω) -h = GT.face_diameter_field(Γ) +#h = GT.face_diameter_field(Γ) GT.vtk_plot(joinpath(outdir,"gamma_ref"),Γref) do plt GT.plot!(plt,g;label="u") GT.plot!(plt,n;label="n") - GT.plot!(plt,q->n2(ϕ(q));label="n2") - GT.plot!(plt,h;label="h") + GT.plot!(plt,q->n2(ϕ(q));label="n2") # TODO + # GT.plot!(plt,h;label="h") TODO GT.plot!(plt;label="u2") do q x = ϕ(q) uref(x) @@ -70,7 +66,8 @@ GT.vtk_plot(joinpath(outdir,"gamma"),Γ) do plt GT.plot!(plt,u;label="u") GT.plot!(plt,n;label="n") GT.plot!(plt,n2;label="n2") - GT.plot!(plt,h;label="h") + # TODO + #GT.plot!(plt,h;label="h") end Λref = GT.skeleton(mesh; @@ -81,38 +78,35 @@ end Λ = GT.physical_domain(Λref) +# TODO n = GT.unit_normal(Λref,Ω) - n2 = GT.unit_normal(Λ,Ω) -h = GT.face_diameter_field(Λ) +#h = GT.face_diameter_field(Λ) -jump(u,q) = u[2](q)-u[1](q) -jump(u,ϕ,q) = u(ϕ[2](q))[2]-u(ϕ[1](q))[1] +jump(u) = u[+] - u[-] GT.vtk_plot(joinpath(outdir,"lambda_ref"),Λref) do plt - GT.plot!(plt,n[1];label="n1") - GT.plot!(plt,n[2];label="n2") + GT.plot!(plt,n[+];label="n1") + GT.plot!(plt,n[-];label="n2") GT.plot!(plt;label="jump_u2") do q - jump(uref∘ϕ,q) + jump((uref∘ϕ)(q)) end - GT.plot!(plt,h;label="h") + #GT.plot!(plt,h;label="h") GT.plot!(plt;label="jump_u") do q - jump(uref,ϕ,q) + jump(uref(ϕ(q))) end end -jump(u) = u[2]-u[1] -jump(u,x) = u(x)[2]-u(x)[1] +jump2(u) = q -> jump(u(q)) + GT.vtk_plot(joinpath(outdir,"lambda"),Λ) do plt - GT.plot!(plt,n2[1];label="n1") - GT.plot!(plt,n2[2];label="n2") + GT.plot!(plt,n2[+];label="n1") + GT.plot!(plt,n2[-];label="n2") GT.plot!(plt;label="jump_u") do q - jump(u,q) - end - GT.plot!(plt;label="jump_u2") do q jump(u(q)) end - GT.plot!(plt,h;label="h") + GT.plot!(plt,jump2(u);label="jump_u2") + #GT.plot!(plt,h;label="h") end # Parallel From 58c58c8cc6b68d2dce254f0b950e0e1a1659fcd6 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Sat, 6 Jul 2024 17:48:21 +0200 Subject: [PATCH 08/20] Deactivating some tests temporarily --- GalerkinToolkitExamples/test/runtests.jl | 2 +- test/runtests.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/GalerkinToolkitExamples/test/runtests.jl b/GalerkinToolkitExamples/test/runtests.jl index 8eefcb73..aecba88b 100644 --- a/GalerkinToolkitExamples/test/runtests.jl +++ b/GalerkinToolkitExamples/test/runtests.jl @@ -3,7 +3,7 @@ module GalerkinToolkitExamplesTests using Test @testset "GalerkinToolkitExamples" begin - @testset "poisson" begin include("poisson_tests.jl") end + #@testset "poisson" begin include("poisson_tests.jl") end end end # module diff --git a/test/runtests.jl b/test/runtests.jl index 0df26b90..24a4bef6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,9 +6,9 @@ using Test @testset "Geometry" begin include("geometry_tests.jl") end @testset "Symbolics" begin include("symbolics_tests.jl") end @testset "Domain" begin include("domain_tests.jl") end - @testset "Integration" begin include("integration_tests.jl") end - @testset "Interpolation" begin include("interpolation_tests.jl") end - @testset "Assembly" begin include("assembly_tests.jl") end + #@testset "Integration" begin include("integration_tests.jl") end + #@testset "Interpolation" begin include("interpolation_tests.jl") end + #@testset "Assembly" begin include("assembly_tests.jl") end end end # module From 9490e68a1e2013e54e3d1566518e3b67d65ab6b0 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Sat, 6 Jul 2024 18:42:43 +0200 Subject: [PATCH 09/20] Starting to add codegen to integration.jl --- src/integration.jl | 137 ++++++++++++++++++++++++-------------- test/integration_tests.jl | 5 +- 2 files changed, 91 insertions(+), 51 deletions(-) diff --git a/src/integration.jl b/src/integration.jl index 8bb93287..d0351a69 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -200,11 +200,15 @@ function coordinates(measure::Measure,domain::ReferenceDomain) prototype = first(GT.coordinates(first(refid_to_quad))) GT.quantity(prototype,domain) do index domface = index.face - face = domface_to_face[domface] - refid = face_to_refid[face] - coords = refid_to_coords[refid] point = index.point - coords[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 end @@ -235,27 +239,34 @@ function weights(measure::Measure,domain::ReferenceDomain) face_to_refid = GT.face_reference_id(mesh,d) domface_to_face = GT.faces(domain) refid_to_quad = reference_quadratures(measure) - refid_to_w = map(GT.weights,refid_to_quad) + 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 - face = domface_to_face[domface] - refid = face_to_refid[face] - w = refid_to_w[refid] point = index.point - w[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 + face = $domface_to_face_sym[$domface] + ws = reference_value($refid_to_ws_sym,$face_to_refid_sym,face) + ws[$point] + end end end +function change_of_measure(J) + sqrt(det(transpose(J)*J)) +end + function weights(measure::Measure,domain::PhysicalDomain) Ω = domain Ωref = GT.reference_domain(Ω) ϕ = GT.domain_map(Ωref,Ω) w = GT.weights(measure,Ωref) x = GT.coordinates(measure,Ωref) - f(J) = sqrt(det(transpose(J)*J)) J = GT.call(ForwardDiff.jacobian,ϕ,x) - dV = GT.call(f,J) + dV = GT.call(change_of_measure,J) GT.call(*,w,dV) end @@ -266,13 +277,17 @@ function num_points(measure::Measure) face_to_refid = GT.face_reference_id(mesh,d) domface_to_face = GT.faces(domain) refid_to_quad = reference_quadratures(measure) - refid_to_w = map(GT.weights,refid_to_quad) + refid_to_ws = map(GT.weights,refid_to_quad) index -> begin domface = index.face - face = domface_to_face[domface] - refid = face_to_refid[face] - w = refid_to_w[refid] - length(w) + 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 + face = $domface_to_face_sym[$domface] + ws = reference_value($refid_to_ws_sym,$face_to_refid_sym,face) + length(ws) + end end end @@ -281,25 +296,21 @@ function integrate(f,measure::Measure) x = GT.coordinates(measure) fx = f(x) contrib = integrate_impl(fx,measure) - integral((domain=>contrib,)) + integral((measure=>contrib,)) end const ∫ = integrate function integrate_impl(fx,measure) w = GT.weights(measure) - p_fx = GT.prototype(fx) - p_w = GT.prototype(w) - t_fx = GT.term(fx) - t_w = GT.term(w) - prototype = p_fx*p_w + p_fx*p_w - num_points = GT.num_points(measure) - contrib = GT.quantity(prototype,GT.domain(measure)) do index - np = num_points(index) - sum(1:np) do point - index2 = replace_point(index,point) - t_fx(index2)*t_w(index2) - end - end + GT.call(*,fx,w) + #p_fx = GT.prototype(fx) + #p_w = GT.prototype(w) + #t_fx = GT.term(fx) + #t_w = GT.term(w) + #prototype = p_fx*p_w + p_fx*p_w + #GT.quantity(prototype,GT.domain(measure)) do index + # GT.call(*,t_fx(index),t_w(index)) + #end end function integrate(f,measure::Measure{<:PMesh}) @@ -310,14 +321,14 @@ function integrate(f,measure::Measure{<:PMesh}) prototype = map(GT.prototype,q) |> PartitionedArrays.getany domain = measure |> GT.domain contrib = GT.quantity(term,prototype,domain) - integral((domain=>contrib,)) + integral((measure=>contrib,)) end function integral(contribs) Integral(contribs) end -# TODO Rename contributions with domain_contribution +# TODO Rename contributions with measure_contribution ? struct Integral{A} contributions::A end @@ -328,39 +339,67 @@ function Base.sum(int::Integral) end function sum_contribution(contrib) - domain, qty = contrib - sum_contribution(domain,qty) + measure, qty = contrib + sum_contribution(measure,qty) end -function sum_contribution(domain,qty) +function sum_contribution(measure,qty) + domain = GT.domain(measure) nfaces = GT.num_faces(domain) facemask = fill(true,nfaces) - sum_contribution_impl(qty,facemask) + sum_contribution_impl(qty,measure,facemask) end -function sum_contribution(domain::AbstractDomain{<:PMesh},qty::AbstractQuantity{<:PMesh}) +function sum_contribution(measure::AbstractDomain{<:PMesh},qty::AbstractQuantity{<:PMesh}) + domain = GT.domain(measure) mesh = domain |> GT.mesh d = GT.num_dims(domain) # TODO allow the user to skip or not to skip ghosts - map(partition(domain),partition(qty),GT.face_partition(mesh,d)) do mydom,myqty,myfaces + map(partition(measure),partition(qty),GT.face_partition(mesh,d)) do mymeasure,myqty,myfaces + mydom = GT.domain(mymeasure) facemask = (part_id(myfaces) .== local_to_owner(myfaces))[GT.faces(mydom)] - sum_contribution_impl(myqty,facemask) + sum_contribution_impl(myqty,mymeasure,facemask) end |> sum end -function sum_contribution_impl(qty,facemask) +function sum_contribution_impl(qty,measure,facemask) # TODO some code duplication with face_contribution_impl - z = zero(GT.prototype(qty)) - nfaces = length(facemask) - term = GT.term(qty) - for face in 1:nfaces - if ! facemask[face] - continue + # 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) |> 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)) + $(s_qty[1]) + $(s_npoints[1]) + nfaces = length(facemask) + for $face in 1:nfaces + if ! facemask[$face] + continue + end + $(s_qty[2]) + npoints = $(s_npoints[2]) + for $point in 1:npoints + s += $(s_qty[3]) + end + end + s end - index = GT.index(;face) - z += term(index) end - z + loop = eval(expr) + storage = GT.storage(index) + z = zero(GT.prototype(qty)) + s = invokelatest(loop,z,facemask,storage) + s end function face_contribution(int::Integral,domain) diff --git a/test/integration_tests.jl b/test/integration_tests.jl index 8d223ec1..408bab0d 100644 --- a/test/integration_tests.jl +++ b/test/integration_tests.jl @@ -53,7 +53,7 @@ int = ∫(dΩref) do q dV = abs(det(J)) u(x)*dV end - +sum(int) @test sum(int) ≈ 8 int = ∫(dΩref) do q @@ -75,7 +75,6 @@ u = GT.analytical_field(x->1,Ω) int = ∫(u,dΩ) @test sum(int) ≈ 4 - D = GT.num_dims(mesh) Γref = GT.boundary(mesh; is_reference_domain=true, @@ -119,6 +118,8 @@ end @test sum(int) ≈ 8 +xxxx + @test sum(GT.face_diameter(Γ)) ≈ 4 h = GT.face_diameter_field(Γ) From 95c8719a7888ae9d02d540226d3c1691d8bdafb0 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Fri, 23 Aug 2024 18:05:25 +0200 Subject: [PATCH 10/20] Fixing integration --- Project.toml | 2 +- src/integration.jl | 61 +++++++++++++++++++++++++++++---------- test/integration_tests.jl | 13 +++++---- test/runtests.jl | 2 +- 4 files changed, 54 insertions(+), 24 deletions(-) diff --git a/Project.toml b/Project.toml index 0ec8a5bb..2a7df090 100644 --- a/Project.toml +++ b/Project.toml @@ -26,7 +26,7 @@ Gmsh = "0.3" IterativeSolvers = "0.9" MPI = "0.20" Metis = "1" -PartitionedArrays = "0.4" +PartitionedArrays = "0.4, 0.5" StaticArrays = "1" WriteVTK = "1" julia = "1" diff --git a/src/integration.jl b/src/integration.jl index d0351a69..f90d4138 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -350,7 +350,7 @@ function sum_contribution(measure,qty) sum_contribution_impl(qty,measure,facemask) end -function sum_contribution(measure::AbstractDomain{<:PMesh},qty::AbstractQuantity{<:PMesh}) +function sum_contribution(measure::Measure{<:PMesh},qty::AbstractQuantity{<:PMesh}) domain = GT.domain(measure) mesh = domain |> GT.mesh d = GT.num_dims(domain) @@ -403,33 +403,62 @@ function sum_contribution_impl(qty,measure,facemask) end function face_contribution(int::Integral,domain) - for (domain2,qty) in int.contributions + for (measure,qty) in int.contributions + domain2 = GT.domain(measure) if domain2 == domain - return face_contribution(domain2,qty) + return face_contribution(measure,qty) end end error("Domain not in integral") end -function face_contribution(domain,qty) +function face_contribution(measure,qty) + domain = GT.domain(measure) nfaces = GT.num_faces(domain) facemask = fill(true,nfaces) - face_contribution_impl(qty,facemask) + face_contribution_impl(qty,measure,facemask) end -function face_contribution_impl(qty,facemask) - z = zero(GT.prototype(qty)) - nfaces = length(facemask) - r = fill(z,nfaces) - term = GT.term(qty) - for face in 1:nfaces - if ! facemask[face] - continue +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 + # 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)) + $(s_qty[1]) + $(s_npoints[1]) + nfaces = length(facemask) + z = zero(eltype(facevals)) + for $face in 1:nfaces + if ! facemask[$face] + continue + end + $(s_qty[2]) + npoints = $(s_npoints[2]) + s = z + for $point in 1:npoints + s += $(s_qty[3]) + end + facevals[$face] = s + end + facevals end - index = GT.index(;face) - r[face] = term(index) end - r + loop! = eval(expr) + storage = GT.storage(index) + z = zero(GT.prototype(qty)) + nfaces = length(facemask) + facevals = fill(z,nfaces) + invokelatest(loop!,facevals,facemask,storage) + facevals end function face_diameter(Ω) diff --git a/test/integration_tests.jl b/test/integration_tests.jl index 408bab0d..38a386c1 100644 --- a/test/integration_tests.jl +++ b/test/integration_tests.jl @@ -118,9 +118,8 @@ end @test sum(int) ≈ 8 -xxxx - -@test sum(GT.face_diameter(Γ)) ≈ 4 +fd = GT.face_diameter(Γ) +@test sum(fd) ≈ 4 h = GT.face_diameter_field(Γ) @@ -133,14 +132,15 @@ dΛref = GT.measure(Λref,degree) ϕ_Λref_Λ = GT.domain_map(Λref,Λ) ϕ_Λref_Ωref = GT.domain_map(Λref,Ωref) -jump(u,ϕ,q) = u(ϕ[2](q))[2]-u(ϕ[1](q))[1] +jump(u,ϕ,q) = u(ϕ[+](q))[+]-u(ϕ[-](q))[-] int = 10*∫(dΛref) do p J = ForwardDiff.jacobian(ϕ_Λref_Λ,p) 3*jump(uref,ϕ_Λref_Ωref,p)*dS(J) end -@test sum(int*1) + 1 ≈ 1 +s = sum(int*1) +@test s + 1 ≈ 1 @test sum(int/1) + 1 ≈ 1 # Parallel @@ -169,7 +169,8 @@ int = ∫(dΩref) do q u(x)*dV end -@test sum(int) ≈ 8 +s = sum(int) +@test s ≈ 8 int = ∫(dΩref) do q J = ForwardDiff.jacobian(ϕ,q) diff --git a/test/runtests.jl b/test/runtests.jl index 24a4bef6..1dbfc103 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ using Test @testset "Geometry" begin include("geometry_tests.jl") end @testset "Symbolics" begin include("symbolics_tests.jl") end @testset "Domain" begin include("domain_tests.jl") end - #@testset "Integration" begin include("integration_tests.jl") end + @testset "Integration" begin include("integration_tests.jl") end #@testset "Interpolation" begin include("interpolation_tests.jl") end #@testset "Assembly" begin include("assembly_tests.jl") end end From 111bec5357c4b3f1e49ca8b1d65fde0c19c8d122 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Sun, 25 Aug 2024 11:06:10 +0200 Subject: [PATCH 11/20] Saving temporary changes [skip ci] --- src/interpolation.jl | 344 +++++++++++++++++++++++++++++++------------ 1 file changed, 252 insertions(+), 92 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index ab67d8da..99347994 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -473,11 +473,19 @@ end function num_face_dofs(a::AbstractSpace,dim,field) face_to_dofs = face_dofs(a) index -> begin + dict = index.dict + field_sym = get!(dict,field,gensym("field")) + face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) + dim_sym = get!(dict,dim,gensym("dim")) face = index.face - if (index.field_per_dim !== nothing && index.field_per_dim[dim] != field) || face == 0 - return 0 + field_per_dim = index.field_per_dim + @assert face !== nothing + @assert field_per_dim !== nothing + @term begin + f = length($face_to_dofs_sym[$face]) + bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face == 0) + ifelse(bool,0,f) end - length(face_to_dofs[face]) end end @@ -489,33 +497,47 @@ end function dof_map(a::AbstractSpace,dim,field) face_to_dofs = face_dofs(a) T = eltype(eltype(face_to_dofs)) + z = T(-1) index -> begin + dict = index.dict + field_sym = get!(dict,field,gensym("field")) + z_sym = get!(dict,z,gensym("z")) + face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) + dim_sym = get!(dict,dim,gensym("dim")) face = index.face - if (index.field_per_dim !== nothing && index.field_per_dim[dim] != field) || face == 0 - return T(-1) + dof = index.dof + dof_per_dim = index.dof_per_dim + field_per_dim = index.field_per_dim + @assert face !== nothing + @assert dof !== nothing + @assert dof_per_dim !== nothing + @assert field_per_dim !== nothing + @term begin + dof = $dof_per_dim[$dim_sym] + f = $face_to_dofs_sym[$face][dof] + bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face == 0) + ifelse(bool,$z_sym,f) end - dof = index.dof_per_dim[dim] - face_to_dofs[face][dof] end end -# TODO I guess this is not needed anymore -function shape_function_mask(f,face_around_per_dim,face_around,dim) - @assert face_around !== nothing - x -> face_around_per_dim[dim] == face_around ? f(x) : zero(f(x)) -end - -function shape_function_mask(f,face_around_per_dim::Nothing,face_around,dim) - f -end - -function shape_function_mask(f,face_around_per_dim::Nothing,face_around::Nothing,dim) - f -end - -function shape_function_mask(f,face_around_per_dim,face_around::Nothing,dim) - f -end +## TODO I guess this is not needed anymore +#function shape_function_mask(f,face_around_per_dim,face_around,dim) +# @assert face_around !== nothing +# x -> face_around_per_dim[dim] == face_around ? f(x) : zero(f(x)) +#end +# +#function shape_function_mask(f,face_around_per_dim::Nothing,face_around,dim) +# f +#end +# +#function shape_function_mask(f,face_around_per_dim::Nothing,face_around::Nothing,dim) +# f +#end +# +#function shape_function_mask(f,face_around_per_dim,face_around::Nothing,dim) +# f +#end function primal_map(a::AbstractSpace) nothing @@ -525,6 +547,7 @@ function dual_map(a::AbstractSpace) nothing end + function shape_functions(a::AbstractSpace,dim) field=1 GT.shape_functions(a,dim,field) @@ -537,16 +560,109 @@ function shape_functions(a::AbstractSpace,dim,field) refid_to_funs = map(GT.shape_functions,refid_to_reffes) domain = GT.domain(a) prototype = first(first(refid_to_funs)) + zero_fun = x -> zero(prototype(x)) + qty = GT.quantity(prototype,domain) do index + dict = index.dict + zero_fun_sym = get!(dict,zero_fun,gensym("zero_fun")) + face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) + refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) + dim_sym = get!(dict,dim,gensym("dim")) + field_sym = get!(dict,field,gensym("field")) + face = index.face + dof = index.dof + dof_per_dim = index.dof_per_dim + field_per_dim = index.field_per_dim + face_around_per_dim = index.face_around_per_dim + face_around = index.face_around + @assert face !== nothing + @assert dof !== nothing + @assert dof_per_dim !== nothing + @assert field_per_dim !== nothing + @assert face_around_per_dim !== nothing + @assert face_around !== nothing + @term begin + refid = $face_to_refid_sym[$face] + dof = $dof_per_dim[$dim_sym] + f = $refid_to_funs_sym[refid][dof] + bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face_around_per_dim[$dim_sym] != $face_around) || ($face == 0) + ifelse(bool,$zero_fun_sym,f) + end + end + if is_reference_domain(GT.domain(a)) + return qty + end + Ω = GT.domain(a) + Ωref = GT.reference_domain(Ω) + ϕ = GT.domain_map(Ωref,Ω) + ϕinv = GT.inverse_map(ϕ) + compose(qty,ϕinv) +end + +function face_function_free_and_dirichlet( + rid_to_fs, + face_to_rid, + face_to_dofs, + free_dofs_to_value, + diri_dofs_to_value, + face, + bool, + z) + if bool + return z + end + fs = reference_value(rid_to_fs,face_to_rid,face) + dofs = face_to_dofs[face] + x -> begin + ndofs = length(dofs) + sum(1:ndofs) do idof + dof = dofs[idof] + coeff = dof > 0 ? free_dofs_to_value[dof] : diri_dofs_to_value[-dof] + f = fs[idof] + f(x)*coeff + end + end +end + +function discrete_field_qty(a::AbstractSpace,free_vals,diri_vals) + @assert primal_map(a) === nothing + face_to_refid = face_reference_id(a) + refid_to_reffes = reference_fes(a) + refid_to_funs = map(GT.shape_functions,refid_to_reffes) + domain = GT.domain(a) + f_prototype = first(first(refid_to_funs)) + z = zero(eltype(free_vals)) + prototype = x-> z*f_prototype(x) + z*f_prototype(x) qty = GT.quantity(prototype,domain) do index + dict = index.dict + zero_fun_sym = get!(dict,zero_fun,gensym("zero_fun")) + face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) + refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) + dim_sym = get!(dict,dim,gensym("dim")) + field_sym = get!(dict,field,gensym("field")) face = index.face - if (index.field_per_dim !== nothing && index.field_per_dim[dim] != field) || face == 0 - return x -> zero(prototype(x)) + dof = index.dof + dof_per_dim = index.dof_per_dim + field_per_dim = index.field_per_dim + face_around_per_dim = index.face_around_per_dim + face_around = index.face_around + @assert face !== nothing + @assert dof !== nothing + @assert dof_per_dim !== nothing + @assert field_per_dim !== nothing + @assert face_around_per_dim !== nothing + @assert face_around !== nothing + @term begin + refid = $face_to_refid_sym[$face] + dof = $dof_per_dim[$dim_sym] + f = $refid_to_funs_sym[refid][dof] + bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face_around_per_dim[$dim_sym] != $face_around) || ($face == 0) + + + + + + ifelse(bool,$zero_fun_sym,f) end - refid = face_to_refid[face] - dof = index.dof_per_dim[dim] - f = refid_to_funs[refid][dof] - g = shape_function_mask(f,index.face_around_per_dim,index.face_around,dim) - shape_function_mask(g,index.field_per_dim,field,dim) end if is_reference_domain(GT.domain(a)) return qty @@ -566,10 +682,21 @@ function dual_basis(a::AbstractSpace,dim) domain = GT.domain(a) prototype = first(first(refid_to_funs)) qty = GT.quantity(prototype,domain) do index + dict = index.dict + face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) + refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) + dim_sym = get!(dict,dim,gensym("dim")) face = index.face - refid = face_to_refid[face] - dof = index.dof_per_dim[dim] - refid_to_funs[refid][dof] + dof = index.dof + dof_per_dim = index.dof_per_dim + @assert face !== nothing + @assert dof !== nothing + @assert dof_per_dim !== nothing + @term begin + refid = $face_to_refid_sym[$face] + dof = $dof_per_dim[$dim_sym] + $refid_to_funs_sym[refid][dof] + end end if is_reference_domain(GT.domain(a)) return qty @@ -702,16 +829,16 @@ function dof_map(a::CartesianProductSpace,dim,field) f = component(a,field) term = dof_map(f,dim,field) index -> begin - dof = term(index) - prolongate_dof(strategy,dof,field) + dof_term = term(index) + dict = index.dict + field_sym = get!(dict,field,gensym("field")) + strategy_sym = get!(dict,strategy,gensym("strategy")) + @term begin + dof = $dof_term + prolongate_dof($strategy_sym,dof,$field_sym) + end + end - #terms = map(s->GT.dof_map(s,dim),a.spaces) - #index -> begin - # field = index.field_per_dim[dim] - # @assert field !== nothing - # index2 = replace_field_per_dim(index,nothing) - # terms[field](index2) - #end end function shape_functions(a::CartesianProductSpace,dim) @@ -731,15 +858,17 @@ function dual_basis(a::CartesianProductSpace) end function discrete_field(space,free_values,dirichlet_values) + qty = discrete_field_qty(space,free_values,dirichlet_values) mesh = space |> GT.mesh - DiscreteField(mesh,space,free_values,dirichlet_values) + DiscreteField(mesh,space,free_values,dirichlet_values,qty) end -struct DiscreteField{A,B,C,D} <: GT.AbstractQuantity{A} +struct DiscreteField{A,B,C,D,E} <: GT.AbstractQuantity{A} mesh::A space::B free_values::C dirichlet_values::D + quantity::E end Base.iterate(m::DiscreteField) = iterate(components(m)) @@ -805,42 +934,47 @@ function domain(u::DiscreteField) end function prototype(u::DiscreteField) - space = GT.space(u) - dim = 2 - basis = GT.shape_functions(space,dim) - f = GT.prototype(basis) - T = eltype(GT.free_values(u)) - z = zero(T) - # TODO use also dirichlet values type? - x -> f(x)*z + f(x)*z + GT.prototype(u.quantity) + #space = GT.space(u) + #dim = 2 + #basis = GT.shape_functions(space,dim) + #f = GT.prototype(basis) + #T = eltype(GT.free_values(u)) + #z = zero(T) + ## TODO use also dirichlet values type? + #x -> f(x)*z + f(x)*z end function term(u::DiscreteField) - space = GT.space(u) - free_vals = GT.free_values(u) - diri_vals = GT.dirichlet_values(u) - dim = 2 - dof_map = GT.dof_map(space,dim) - num_face_dofs = GT.num_face_dofs(space,dim) - basis = GT.shape_functions(space,dim) - basis_term = GT.term(basis) - index -> begin - x -> begin - index2 = GT.index( - face=index.face, - local_face=index.local_face, - face_around=index.face_around) - ndofs = num_face_dofs(index2) - sum(1:ndofs) do dof - dof_per_dim = (nothing,dof) - index3 = replace_dof_per_dim(index2,dof_per_dim) - dof = dof_map(index3) - f = basis_term(index3) - coeff = dof > 0 ? free_vals[dof] : diri_vals[-dof] - f(x)*coeff - end - end - end + GT.term(u.quantity) + #space = GT.space(u) + #free_vals = GT.free_values(u) + #diri_vals = GT.dirichlet_values(u) + #space = GT.space(u) + #free_vals = GT.free_values(u) + #diri_vals = GT.dirichlet_values(u) + #dim = 2 + #dof_map = GT.dof_map(space,dim) + #num_face_dofs = GT.num_face_dofs(space,dim) + #basis = GT.shape_functions(space,dim) + #basis_term = GT.term(basis) + #index -> begin + # x -> begin + # index2 = GT.index( + # face=index.face, + # local_face=index.local_face, + # face_around=index.face_around) + # ndofs = num_face_dofs(index2) + # sum(1:ndofs) do dof + # dof_per_dim = (nothing,dof) + # index3 = replace_dof_per_dim(index2,dof_per_dim) + # dof = dof_map(index3) + # f = basis_term(index3) + # coeff = dof > 0 ? free_vals[dof] : diri_vals[-dof] + # f(x)*coeff + # end + # end + #end end # TODO @@ -888,24 +1022,50 @@ function interpolate_impl!(f,u,free_or_diri;location=1) domain = GT.domain(space) num_face_dofs = GT.num_face_dofs(space,dim) dirichlet_dof_location = GT.dirichlet_dof_location(space) - for face in 1:GT.num_faces(domain) - index = GT.index(;face) - for dof in 1:num_face_dofs(index) - index2 = replace_dof_per_dim(index,(dof,)) - v = vals_term(index2) - dof = dof_map(index2) - if dof > 0 - if free_or_diri != DIRICHLET - free_vals[dof] = v - end - else - diri_dof = -dof - if free_or_diri != FREE && dirichlet_dof_location[diri_dof] == location - diri_vals[diri_dof] = v + face = :face + dof = :dof + index = GT.index(;face,dof_per_dim(dof,)) + expr_qty = vals_term(index) |> simplify + expr_ndofs = num_face_dofs(index) |> simplify + expr_map = dof_map(index) |> simplify + # TODO merge statements + s_qty = GT.topological_sort(expr_qty,(face,dof)) + s_map = GT.topological_sort(expr_map,(face,dof)) + s_ndofs = GT.topological_sort(expr_ndofs,(face,)) + expr = quote + function loop!(args,storage) + (;free_vals,diri_vals,nfaces,dirichlet_dof_location,location) = args + $(unpack_storage(index.dict,:storage)) + $(s_qty[1]) + $(s_map[1]) + $(s_ndofs[1]) + for face in 1:nfaces + $(s_qty[2]) + $(s_map[2]) + ndofs = $(s_ndofs[2]) + for idof in 1:ndofs + v = $(s_qty[3]) + dof = $(s_map[3]) + if dof > 0 + if free_or_diri != DIRICHLET + free_vals[dof] = v + end + else + diri_dof = -dof + if free_or_diri != FREE && dirichlet_dof_location[diri_dof] == location + diri_vals[diri_dof] = v + end + end end end end end + loop! = eval(expr) + storage = GT.storage(index) + nfaces = GT.num_faces(domain) + args = (;free_vals,diri_vals,nfaces,dirichlet_dof_location,location) + invokelatest(loop!,args,storage) + u end function interpolate_impl!(f::PiecewiseField,u,free_or_diri) From 2b21cb5eeb1a508ab7d370a382773339a723f310 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 27 Aug 2024 17:58:51 +0200 Subject: [PATCH 12/20] Fixing interpolation --- src/domain.jl | 1 + src/interpolation.jl | 332 ++++++++++++++++++++++--------------------- test/runtests.jl | 2 +- 3 files changed, 173 insertions(+), 162 deletions(-) diff --git a/src/domain.jl b/src/domain.jl index 37df0e78..b15a3028 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -650,6 +650,7 @@ end # end #end +# TODO a better name function face_function(rid_to_fs,face_to_rid,face_to_dofs,dofs_to_value,face) fs = reference_value(rid_to_fs,face_to_rid,face) dofs = face_to_dofs[face] diff --git a/src/interpolation.jl b/src/interpolation.jl index 99347994..4949e0f3 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -455,6 +455,11 @@ function face_dofs(space::AbstractSpace) state.cell_to_dofs # TODO rename face_dofs ? end +function face_dofs(space::AbstractSpace,field) + @assert field == 1 + face_dofs(space) +end + function free_and_dirichlet_dofs(V::AbstractSpace) state = generate_dof_ids(V) state.free_and_dirichlet_dofs @@ -465,61 +470,61 @@ function dirichlet_dof_location(V::AbstractSpace) state.dirichlet_dof_location end -function num_face_dofs(a::AbstractSpace,dim) - field = 1 - num_face_dofs(a,dim,field) -end - -function num_face_dofs(a::AbstractSpace,dim,field) - face_to_dofs = face_dofs(a) - index -> begin - dict = index.dict - field_sym = get!(dict,field,gensym("field")) - face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) - dim_sym = get!(dict,dim,gensym("dim")) - face = index.face - field_per_dim = index.field_per_dim - @assert face !== nothing - @assert field_per_dim !== nothing - @term begin - f = length($face_to_dofs_sym[$face]) - bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face == 0) - ifelse(bool,0,f) - end - end -end - -function dof_map(a::AbstractSpace,dim) - field = 1 - dof_map(a,dim,field) -end - -function dof_map(a::AbstractSpace,dim,field) - face_to_dofs = face_dofs(a) - T = eltype(eltype(face_to_dofs)) - z = T(-1) - index -> begin - dict = index.dict - field_sym = get!(dict,field,gensym("field")) - z_sym = get!(dict,z,gensym("z")) - face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) - dim_sym = get!(dict,dim,gensym("dim")) - face = index.face - dof = index.dof - dof_per_dim = index.dof_per_dim - field_per_dim = index.field_per_dim - @assert face !== nothing - @assert dof !== nothing - @assert dof_per_dim !== nothing - @assert field_per_dim !== nothing - @term begin - dof = $dof_per_dim[$dim_sym] - f = $face_to_dofs_sym[$face][dof] - bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face == 0) - ifelse(bool,$z_sym,f) - end - end -end +#function num_face_dofs(a::AbstractSpace,dim) +# field = 1 +# num_face_dofs(a,dim,field) +#end +# +#function num_face_dofs(a::AbstractSpace,dim,field) +# face_to_dofs = face_dofs(a) +# index -> begin +# dict = index.dict +# field_sym = get!(dict,field,gensym("field")) +# face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) +# dim_sym = get!(dict,dim,gensym("dim")) +# face = index.face +# field_per_dim = index.field_per_dim +# @assert face !== nothing +# @assert field_per_dim !== nothing +# @term begin +# f = length($face_to_dofs_sym[$face]) +# bool = ($field_per_dim[$dim_sym] != $field_sym) | ($face == 0) +# ifelse(bool,0,f) +# end +# end +#end +# +#function dof_map(a::AbstractSpace,dim) +# field = 1 +# dof_map(a,dim,field) +#end +# +#function dof_map(a::AbstractSpace,dim,field) +# face_to_dofs = face_dofs(a) +# T = eltype(eltype(face_to_dofs)) +# z = T(-1) +# index -> begin +# dict = index.dict +# field_sym = get!(dict,field,gensym("field")) +# z_sym = get!(dict,z,gensym("z")) +# face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) +# dim_sym = get!(dict,dim,gensym("dim")) +# face = index.face +# dof = index.dof +# dof_per_dim = index.dof_per_dim +# field_per_dim = index.field_per_dim +# @assert face !== nothing +# @assert dof !== nothing +# @assert dof_per_dim !== nothing +# @assert field_per_dim !== nothing +# @term begin +# dof = $dof_per_dim[$dim_sym] +# f = $face_to_dofs_sym[$face][dof] +# bool = ($field_per_dim[$dim_sym] != $field_sym) | ($face == 0) +# ifelse(bool,$z_sym,f) +# end +# end +#end ## TODO I guess this is not needed anymore #function shape_function_mask(f,face_around_per_dim,face_around,dim) @@ -547,45 +552,49 @@ function dual_map(a::AbstractSpace) nothing end +function shape_functions(a::AbstractSpace,dim,field) + @assert field == 1 + GT.shape_functions(a,dim) +end -function shape_functions(a::AbstractSpace,dim) - field=1 - GT.shape_functions(a,dim,field) +function mask_function(f,bool) + x->bool*f(x) end -function shape_functions(a::AbstractSpace,dim,field) +# TODO a better name +function face_shape_function(rid_to_fs,face_to_rid,face,dof) + fs = reference_value(rid_to_fs,face_to_rid,face) + fs[dof] +end + +function shape_functions(a::AbstractSpace,dim) @assert primal_map(a) === nothing face_to_refid = face_reference_id(a) refid_to_reffes = reference_fes(a) refid_to_funs = map(GT.shape_functions,refid_to_reffes) domain = GT.domain(a) prototype = first(first(refid_to_funs)) - zero_fun = x -> zero(prototype(x)) qty = GT.quantity(prototype,domain) do index dict = index.dict - zero_fun_sym = get!(dict,zero_fun,gensym("zero_fun")) face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) - dim_sym = get!(dict,dim,gensym("dim")) - field_sym = get!(dict,field,gensym("field")) face = index.face - dof = index.dof dof_per_dim = index.dof_per_dim - field_per_dim = index.field_per_dim face_around_per_dim = index.face_around_per_dim face_around = index.face_around @assert face !== nothing - @assert dof !== nothing @assert dof_per_dim !== nothing - @assert field_per_dim !== nothing - @assert face_around_per_dim !== nothing - @assert face_around !== nothing - @term begin - refid = $face_to_refid_sym[$face] - dof = $dof_per_dim[$dim_sym] - f = $refid_to_funs_sym[refid][dof] - bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face_around_per_dim[$dim_sym] != $face_around) || ($face == 0) - ifelse(bool,$zero_fun_sym,f) + dof = dof_per_dim[dim] + if face_around_per_dim !== nothing + @assert face_around === nothing + @term face_shape_function($refid_to_funs_sym,$face_to_refid_sym,$face,$dof) + else + @assert face_around !== nothing + @term begin + fs = face_shape_function($refid_to_funs_sym,$face_to_refid_sym,$face,$dof) + bool = $face_around == $(face_around_per_dim[dim]) + mask_function(fs,bool) + end end end if is_reference_domain(GT.domain(a)) @@ -598,18 +607,14 @@ function shape_functions(a::AbstractSpace,dim,field) compose(qty,ϕinv) end +# TODO a better name function face_function_free_and_dirichlet( rid_to_fs, face_to_rid, face_to_dofs, free_dofs_to_value, diri_dofs_to_value, - face, - bool, - z) - if bool - return z - end + face) fs = reference_value(rid_to_fs,face_to_rid,face) dofs = face_to_dofs[face] x -> begin @@ -632,36 +637,24 @@ function discrete_field_qty(a::AbstractSpace,free_vals,diri_vals) f_prototype = first(first(refid_to_funs)) z = zero(eltype(free_vals)) prototype = x-> z*f_prototype(x) + z*f_prototype(x) + face_to_dofs = GT.face_dofs(a) qty = GT.quantity(prototype,domain) do index dict = index.dict - zero_fun_sym = get!(dict,zero_fun,gensym("zero_fun")) face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) - dim_sym = get!(dict,dim,gensym("dim")) - field_sym = get!(dict,field,gensym("field")) + face_to_dofs_sym = get!(dict,face_to_dofs,gensym("face_to_dofs")) + free_vals_sym = get!(dict,free_vals,gensym("free_vals")) + diri_vals_sym = get!(dict,diri_vals,gensym("diri_vals")) face = index.face - dof = index.dof - dof_per_dim = index.dof_per_dim - field_per_dim = index.field_per_dim - face_around_per_dim = index.face_around_per_dim - face_around = index.face_around @assert face !== nothing - @assert dof !== nothing - @assert dof_per_dim !== nothing - @assert field_per_dim !== nothing - @assert face_around_per_dim !== nothing - @assert face_around !== nothing @term begin - refid = $face_to_refid_sym[$face] - dof = $dof_per_dim[$dim_sym] - f = $refid_to_funs_sym[refid][dof] - bool = ($field_per_dim[$dim_sym] != $field_sym) || ($face_around_per_dim[$dim_sym] != $face_around) || ($face == 0) - - - - - - ifelse(bool,$zero_fun_sym,f) + face_function_free_and_dirichlet( + $refid_to_funs_sym, + $face_to_refid_sym, + $face_to_dofs_sym, + $free_vals_sym, + $diri_vals_sym, + $face) end end if is_reference_domain(GT.domain(a)) @@ -685,17 +678,14 @@ function dual_basis(a::AbstractSpace,dim) dict = index.dict face_to_refid_sym = get!(dict,face_to_refid,gensym("face_to_refid")) refid_to_funs_sym = get!(dict,refid_to_funs,gensym("refid_to_funs")) - dim_sym = get!(dict,dim,gensym("dim")) face = index.face - dof = index.dof dof_per_dim = index.dof_per_dim @assert face !== nothing - @assert dof !== nothing @assert dof_per_dim !== nothing + dof = dof_per_dim[dim] @term begin refid = $face_to_refid_sym[$face] - dof = $dof_per_dim[$dim_sym] - $refid_to_funs_sym[refid][dof] + $refid_to_funs_sym[refid][$dof] end end if is_reference_domain(GT.domain(a)) @@ -804,41 +794,51 @@ function domain(a::CartesianProductSpace,field) GT.domain(GT.component(a,field)) end -function num_face_dofs(a::CartesianProductSpace,dim) - error("Casenot implemented (perhaps not needed in practice)") -end +#function num_face_dofs(a::CartesianProductSpace,dim) +# error("Casenot implemented (perhaps not needed in practice)") +#end +# +#function num_face_dofs(a::CartesianProductSpace,dim,field) +# f = component(a,field) +# num_face_dofs(f,dim,field) +# #terms = map(s->GT.num_face_dofs(s,dim),a.spaces) +# #index -> begin +# # field = index.field_per_dim[dim] +# # @assert field !== nothing +# # index2 = replace_field_per_dim(index,nothing) +# # terms[field](index2) +# #end +#end -function num_face_dofs(a::CartesianProductSpace,dim,field) - f = component(a,field) - num_face_dofs(f,dim,field) - #terms = map(s->GT.num_face_dofs(s,dim),a.spaces) - #index -> begin - # field = index.field_per_dim[dim] - # @assert field !== nothing - # index2 = replace_field_per_dim(index,nothing) - # terms[field](index2) - #end -end +## TODO face_dofs vs dof_map very confusing. Remove dof_map? +#function dof_map(a::CartesianProductSpace,dim) +# error("Casenot implemented (perhaps not needed in practice)") +#end +# +#function dof_map(a::CartesianProductSpace,dim,field) +# strategy = GT.free_values_strategy(a) +# f = component(a,field) +# term = dof_map(f,dim,field) +# index -> begin +# dof_term = term(index) +# dict = index.dict +# field_sym = get!(dict,field,gensym("field")) +# strategy_sym = get!(dict,strategy,gensym("strategy")) +# @term begin +# dof = $dof_term +# prolongate_dof($strategy_sym,dof,$field_sym) +# end +# +# end +#end +# -function dof_map(a::CartesianProductSpace,dim) - error("Casenot implemented (perhaps not needed in practice)") +function face_dofs(a::CartesianProductSpace) + error("Not implemented, not needed in practice") end -function dof_map(a::CartesianProductSpace,dim,field) - strategy = GT.free_values_strategy(a) - f = component(a,field) - term = dof_map(f,dim,field) - index -> begin - dof_term = term(index) - dict = index.dict - field_sym = get!(dict,field,gensym("field")) - strategy_sym = get!(dict,strategy,gensym("strategy")) - @term begin - dof = $dof_term - prolongate_dof($strategy_sym,dof,$field_sym) - end - - end +function face_dofs(a::CartesianProductSpace,field) + face_dofs(a.spaces[field]) end function shape_functions(a::CartesianProductSpace,dim) @@ -850,13 +850,31 @@ end function shape_functions(a::CartesianProductSpace,dim,field) f = component(a,field) - shape_functions(f,dim,field) + qty = shape_functions(f,dim,field) + t = GT.term(qty) + GT.quantity(GT.prototype(qty),GT.domain(qty)) do index + field_per_dim = index.field_per_dim + @assert field_per_dim !== nothing + @term begin + bool = $field == $(field_per_dim[field]) + f = $(t(index)) + mask_function(f,bool) + end + end +end + +function discrete_field_qty(a::CartesianProductSpace,free_vals,diri_vals) + nothing end function dual_basis(a::CartesianProductSpace) error("Not implemented yet. Not needed in practice.") end +function dual_basis(a::CartesianProductSpace,field) + error("Not implemented yet. Not needed in practice.") +end + function discrete_field(space,free_values,dirichlet_values) qty = discrete_field_qty(space,free_values,dirichlet_values) mesh = space |> GT.mesh @@ -1016,42 +1034,34 @@ function interpolate_impl!(f,u,free_or_diri;location=1) space = GT.space(u) dim = 1 sigma = GT.dual_basis(space,dim) - dof_map = GT.dof_map(space,dim) + face_to_dofs = GT.face_dofs(space) vals = sigma(f) vals_term = GT.term(vals) domain = GT.domain(space) - num_face_dofs = GT.num_face_dofs(space,dim) dirichlet_dof_location = GT.dirichlet_dof_location(space) face = :face dof = :dof - index = GT.index(;face,dof_per_dim(dof,)) + index = GT.index(;face,dof_per_dim=(dof,)) expr_qty = vals_term(index) |> simplify - expr_ndofs = num_face_dofs(index) |> simplify - expr_map = dof_map(index) |> simplify - # TODO merge statements s_qty = GT.topological_sort(expr_qty,(face,dof)) - s_map = GT.topological_sort(expr_map,(face,dof)) - s_ndofs = GT.topological_sort(expr_ndofs,(face,)) expr = quote function loop!(args,storage) - (;free_vals,diri_vals,nfaces,dirichlet_dof_location,location) = args + (;face_to_dofs,free_vals,diri_vals,nfaces,dirichlet_dof_location,location,free_or_diri) = args $(unpack_storage(index.dict,:storage)) $(s_qty[1]) - $(s_map[1]) - $(s_ndofs[1]) - for face in 1:nfaces + for $face in 1:nfaces $(s_qty[2]) - $(s_map[2]) - ndofs = $(s_ndofs[2]) - for idof in 1:ndofs + dofs = face_to_dofs[$face] + ndofs = length(dofs) + for $dof in 1:ndofs v = $(s_qty[3]) - dof = $(s_map[3]) - if dof > 0 + gdof = dofs[$dof] + if gdof > 0 if free_or_diri != DIRICHLET - free_vals[dof] = v + free_vals[gdof] = v end else - diri_dof = -dof + diri_dof = -gdof if free_or_diri != FREE && dirichlet_dof_location[diri_dof] == location diri_vals[diri_dof] = v end @@ -1063,7 +1073,7 @@ function interpolate_impl!(f,u,free_or_diri;location=1) loop! = eval(expr) storage = GT.storage(index) nfaces = GT.num_faces(domain) - args = (;free_vals,diri_vals,nfaces,dirichlet_dof_location,location) + args = (;face_to_dofs,free_vals,diri_vals,nfaces,dirichlet_dof_location,location,free_or_diri) invokelatest(loop!,args,storage) u end diff --git a/test/runtests.jl b/test/runtests.jl index 1dbfc103..efe3ebf8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,7 @@ using Test @testset "Symbolics" begin include("symbolics_tests.jl") end @testset "Domain" begin include("domain_tests.jl") end @testset "Integration" begin include("integration_tests.jl") end - #@testset "Interpolation" begin include("interpolation_tests.jl") end + @testset "Interpolation" begin include("interpolation_tests.jl") end #@testset "Assembly" begin include("assembly_tests.jl") end end From 0c6b50b48fc0aa46079adf2a4e391955b051bd07 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 28 Aug 2024 14:12:27 +0200 Subject: [PATCH 13/20] Assembly working with code generation --- src/assembly.jl | 225 +++++++++++++++++++++++++++-------------- src/interpolation.jl | 12 +-- test/assembly_tests.jl | 2 +- test/runtests.jl | 2 +- 4 files changed, 154 insertions(+), 87 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 4d8ca001..33197c58 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -34,8 +34,9 @@ function assemble_vector_count(state) contributions = GT.contributions(integral) num_fields = GT.num_fields(space) fields = ntuple(identity,num_fields) - glues = map(contributions) do domain_and_contribution - domain, = domain_and_contribution + glues = map(contributions) do measure_and_contribution + measure, = measure_and_contribution + domain = GT.domain(measure) map(fields) do field GT.domain_glue(domain,GT.domain(space,field);strict=false) end @@ -45,14 +46,12 @@ function assemble_vector_count(state) end function loop(glue,field) sface_to_faces, = target_face(glue) - dim = 1 - num_face_dofs = GT.num_face_dofs(space,dim,field) - field_per_dim = (field,) - for sface in 1:GT.num_faces(GT.domain(glue)) + face_to_dofs = face_dofs(space,field) + nsfaces = length(sface_to_faces) + for sface in 1:nsfaces faces = sface_to_faces[sface] for face in faces - index = GT.index(;face,field_per_dim) - ndofs = num_face_dofs(index) + ndofs = length(face_to_dofs[face]) n += ndofs end end @@ -76,37 +75,69 @@ function assemble_vector_fill(state) contributions = GT.contributions(integral) num_fields = GT.num_fields(space) n = 0 - function loop(glue::Nothing,field,contribution) + function loop(glue::Nothing,field,measure,contribution) end - function loop(glue,field,contribution) + function loop(glue,field,measure,contribution) sface_to_faces,_,sface_to_faces_around = target_face(glue) - dim =1 - dof_map = GT.dof_map(space,dim,field) - num_face_dofs = GT.num_face_dofs(space,dim,field) + face_to_dofs = face_dofs(space,field) + term_qty = GT.term(contribution) + term_npoints = GT.num_points(measure) + sface = :sface + point = :point + idof = :idof + face_around = :face_around field_per_dim = (field,) - for sface in 1:GT.num_faces(GT.domain(glue)) - faces = sface_to_faces[sface] - faces_around = sface_to_faces_around[sface] - for (face_around,face) in zip(faces_around,faces) - face_around_per_dim = (face_around,) - index = GT.index(;face,field_per_dim) - ndofs = num_face_dofs(index) - for dof in 1:ndofs - n += 1 - dof_per_dim = (dof,) - index2 = GT.index(;face,dof_per_dim,field_per_dim,face_around_per_dim) - I[n] = dof_map(index2) - index3 = GT.replace_face(index2,sface) - vn = GT.term(contribution)(index3) - V[n] = vn + dof_per_dim = (idof,) + face_around_per_dim = (face_around,) + index = GT.index(;face=sface,point,field_per_dim,dof_per_dim,face_around_per_dim) + expr_qty = term_qty(index) |> simplify + expr_npoints = term_npoints(index) |> simplify + # TODO merge statements + s_qty = GT.topological_sort(expr_qty,(sface,face_around,idof,point)) + s_npoints = GT.topological_sort(expr_npoints,(sface,)) + expr = quote + function loop!(n,args,storage) + (;sface_to_faces,sface_to_faces_around,face_to_dofs,I,V) = args + $(unpack_storage(index.dict,:storage)) + $(s_qty[1]) + $(s_npoints[1]) + nsfaces = length(sface_to_faces) + z = zero(eltype(V)) + for $sface in 1:nsfaces + $(s_qty[2]) + npoints = $(s_npoints[2]) + faces = sface_to_faces[sface] + faces_around = sface_to_faces_around[sface] + for ($face_around,face) in zip(faces_around,faces) + $(s_qty[3]) + dofs = face_to_dofs[face] + ndofs = length(dofs) + for $idof in 1:ndofs + $(s_qty[4]) + v = z + for $point in 1:npoints + v += $(s_qty[5]) + end + n += 1 + dof = dofs[$idof] + I[n] = dof + V[n] = v + end + end end + n end end + loop! = eval(expr) + args = (;sface_to_faces,sface_to_faces_around,face_to_dofs,I,V) + storage = GT.storage(index) + n = invokelatest(loop!,n,args,storage) + nothing end - for (domain_and_contribution,field_to_glue) in zip(contributions,glues) - domain, contribution = domain_and_contribution + for (measure_and_contribution,field_to_glue) in zip(contributions,glues) + measure, contribution = measure_and_contribution map(fields,field_to_glue) do field,glue - loop(glue,field,contribution) + loop(glue,field,measure,contribution) end end state @@ -150,37 +181,36 @@ function assemble_matrix_count(state) num_fields_trial = GT.num_fields(trial_space) fields_test = ntuple(identity,num_fields_test) fields_trial = ntuple(identity,num_fields_trial) - glues_test = map(contributions) do domain_and_contribution - domain, = domain_and_contribution + glues_test = map(contributions) do measure_and_contribution + measure, = measure_and_contribution + domain = GT.domain(measure) map(fields_test) do field GT.domain_glue(domain,GT.domain(test_space,field),strict=false) end end - glues_trial = map(contributions) do domain_and_contribution - domain, = domain_and_contribution + glues_trial = map(contributions) do measure_and_contribution + measure, = measure_and_contribution + domain = GT.domain(measure) map(fields_trial) do field GT.domain_glue(domain,GT.domain(trial_space,field),strict=false) end end n = 0 - # TODO some cross terms missing - # TODO a lot of code duplication function loop(glue_test,glue_trial,field_per_dim) sface_to_faces_test, = target_face(glue_test) sface_to_faces_trial, = target_face(glue_trial) - test_dim = 1 - trial_dim = 2 - num_face_dofs_test = GT.num_face_dofs(test_space,test_dim,field_per_dim[1]) - num_face_dofs_trial = GT.num_face_dofs(trial_space,trial_dim,field_per_dim[2]) - for sface in 1:GT.num_faces(GT.domain(glue_test)) + face_to_dofs_test = face_dofs(test_space,field_per_dim[1]) + face_to_dofs_trial = face_dofs(trial_space,field_per_dim[2]) + nsfaces = length(sface_to_faces_test) + for sface in 1:nsfaces faces_test = sface_to_faces_test[sface] faces_trial = sface_to_faces_trial[sface] for face_test in faces_test - index_test = GT.index(;face=face_test,field_per_dim) - ndofs_test = num_face_dofs_test(index_test) + dofs_test = face_to_dofs_test[face_test] + ndofs_test = length(dofs_test) for face_trial in faces_trial - index_trial = GT.index(;face=face_trial,field_per_dim) - ndofs_trial = num_face_dofs_trial(index_trial) + dofs_trial = face_to_dofs_trial[face_trial] + ndofs_trial = length(dofs_trial) n += ndofs_test*ndofs_trial end end @@ -217,45 +247,84 @@ function assemble_matrix_fill(state) num_fields_test = GT.num_fields(test_space) num_fields_trial = GT.num_fields(trial_space) n = 0 - function loop(glue_test,glue_trial,field_per_dim,contribution) + function loop(glue_test,glue_trial,field_per_dim,measure,contribution) sface_to_faces_test, _,sface_to_faces_around_test = target_face(glue_test) sface_to_faces_trial, _,sface_to_faces_around_trial = target_face(glue_trial) - test_dim = 1 - trial_dim = 2 - num_face_dofs_test = GT.num_face_dofs(test_space,test_dim,field_per_dim[1]) - num_face_dofs_trial = GT.num_face_dofs(trial_space,trial_dim,field_per_dim[2]) - dof_map_test = GT.dof_map(test_space,test_dim,field_per_dim[1]) - dof_map_trial = GT.dof_map(trial_space,trial_dim,field_per_dim[2]) - for sface in 1:GT.num_faces(GT.domain(glue_test)) - faces_test = sface_to_faces_test[sface] - faces_trial = sface_to_faces_trial[sface] - faces_around_test = sface_to_faces_around_test[sface] - faces_around_trial = sface_to_faces_around_trial[sface] - for (face_around_test,face_test) in zip(faces_around_test,faces_test) - index_test = GT.index(;face=face_test,field_per_dim) - ndofs_test = num_face_dofs_test(index_test) - for (face_around_trial,face_trial) in zip(faces_around_trial,faces_trial) - face_around_per_dim = (face_around_test,face_around_trial) - index_trial = GT.index(;face=face_trial,field_per_dim) - ndofs_trial = num_face_dofs_trial(index_trial) - for dof_test in 1:ndofs_test - for dof_trial in 1:ndofs_trial - n += 1 - dof_per_dim = (dof_test,dof_trial) - index2_test = GT.index(;face=face_test,dof_per_dim,field_per_dim,face_around_per_dim) - index2_trial = GT.index(;face=face_trial,dof_per_dim,field_per_dim,face_around_per_dim) - I[n] = dof_map_test(index2_test) - J[n] = dof_map_trial(index2_trial) - index3 = GT.replace_face(index2_test,sface) - V[n] = GT.term(contribution)(index3) + face_to_dofs_test = face_dofs(test_space,field_per_dim[1]) + face_to_dofs_trial = face_dofs(trial_space,field_per_dim[2]) + term_qty = GT.term(contribution) + term_npoints = GT.num_points(measure) + sface = :sface + point = :point + idof_test = :idof_test + idof_trial = :idof_trial + face_around_test = :face_around_test + face_around_trial = :face_around_trial + dof_per_dim = (idof_test,idof_trial) + face_around_per_dim = (face_around_test,face_around_trial) + index = GT.index(;face=sface,point,field_per_dim,dof_per_dim,face_around_per_dim) + expr_qty = term_qty(index) |> simplify + expr_npoints = term_npoints(index) |> simplify + # TODO merge statements + s_qty = GT.topological_sort(expr_qty,(sface,face_around_test,face_around_trial,idof_test,idof_trial,point)) + s_npoints = GT.topological_sort(expr_npoints,(sface,)) + expr = quote + function loop!(n,args,storage) + (;sface_to_faces_test,sface_to_faces_around_test, + sface_to_faces_trial,sface_to_faces_around_trial, + face_to_dofs_test,face_to_dofs_trial,I,J,V) = args + $(unpack_storage(index.dict,:storage)) + $(s_qty[1]) + $(s_npoints[1]) + nsfaces = length(sface_to_faces_test) + z = zero(eltype(V)) + for $sface in 1:nsfaces + $(s_qty[2]) + npoints = $(s_npoints[2]) + faces_test = sface_to_faces_test[sface] + faces_trial = sface_to_faces_trial[sface] + faces_around_test = sface_to_faces_around_test[sface] + faces_around_trial = sface_to_faces_around_trial[sface] + for ($face_around_test,face_test) in zip(faces_around_test,faces_test) + $(s_qty[3]) + dofs_test = face_to_dofs_test[face_test] + ndofs_test = length(dofs_test) + for ($face_around_trial,face_trial) in zip(faces_around_trial,faces_trial) + $(s_qty[4]) + dofs_trial = face_to_dofs_trial[face_trial] + ndofs_trial = length(dofs_trial) + for $idof_test in 1:ndofs_test + $(s_qty[5]) + dof_test = dofs_test[$idof_test] + for $idof_trial in 1:ndofs_trial + $(s_qty[6]) + dof_trial = dofs_trial[$idof_trial] + v = z + for $point in 1:npoints + v += $(s_qty[7]) + end + n += 1 + I[n] = dof_test + J[n] = dof_trial + V[n] = v + end + end end end end + n end end + loop! = eval(expr) + storage = GT.storage(index) + args = (;sface_to_faces_test,sface_to_faces_around_test, + sface_to_faces_trial,sface_to_faces_around_trial, + face_to_dofs_test,face_to_dofs_trial,I,J,V) + n = invokelatest(loop!,n,args,storage) + nothing end - for (domain_and_contribution,field_to_glue_test,field_to_glue_trial) in zip(contributions,glues_test,glues_trial) - domain, contribution = domain_and_contribution + for (measure_and_contribution,field_to_glue_test,field_to_glue_trial) in zip(contributions,glues_test,glues_trial) + measure, contribution = measure_and_contribution for field_test in fields_test glue_test = field_to_glue_test[field_test] for field_trial in fields_trial @@ -264,7 +333,7 @@ function assemble_matrix_fill(state) continue end field_per_dim = (field_test,field_trial) - loop(glue_test,glue_trial,field_per_dim,contribution) + loop(glue_test,glue_trial,field_per_dim,measure,contribution) end end end diff --git a/src/interpolation.jl b/src/interpolation.jl index 4949e0f3..e6a3e06f 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -585,11 +585,9 @@ function shape_functions(a::AbstractSpace,dim) @assert face !== nothing @assert dof_per_dim !== nothing dof = dof_per_dim[dim] - if face_around_per_dim !== nothing - @assert face_around === nothing + if face_around_per_dim === nothing || face_around == nothing @term face_shape_function($refid_to_funs_sym,$face_to_refid_sym,$face,$dof) else - @assert face_around !== nothing @term begin fs = face_shape_function($refid_to_funs_sym,$face_to_refid_sym,$face,$dof) bool = $face_around == $(face_around_per_dim[dim]) @@ -850,15 +848,15 @@ end function shape_functions(a::CartesianProductSpace,dim,field) f = component(a,field) - qty = shape_functions(f,dim,field) + qty = shape_functions(f,dim) t = GT.term(qty) GT.quantity(GT.prototype(qty),GT.domain(qty)) do index + # TODO field_per_dim[dim] has possibly already a numeric value field_per_dim = index.field_per_dim @assert field_per_dim !== nothing @term begin - bool = $field == $(field_per_dim[field]) - f = $(t(index)) - mask_function(f,bool) + bool = $field == $(field_per_dim[dim]) + mask_function($(t(index)),bool) end end end diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index 749ee6e0..e9b56e21 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -55,7 +55,7 @@ function dS(J) sqrt(det(Jt*J)) end -jump(u,ϕ,q) = u(ϕ[2](q))[2]-u(ϕ[1](q))[1] +jump(u,ϕ,q) = u(ϕ[+](q))[+]-u(ϕ[-](q))[-] function l(v) ∫(dΩref) do q diff --git a/test/runtests.jl b/test/runtests.jl index efe3ebf8..0df26b90 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,7 @@ using Test @testset "Domain" begin include("domain_tests.jl") end @testset "Integration" begin include("integration_tests.jl") end @testset "Interpolation" begin include("interpolation_tests.jl") end - #@testset "Assembly" begin include("assembly_tests.jl") end + @testset "Assembly" begin include("assembly_tests.jl") end end end # module From 82b85370ec1bb99351c26578929fc3e92fd2996d Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 28 Aug 2024 17:19:43 +0200 Subject: [PATCH 14/20] Added vector assembly strategy --- Project.toml | 2 + src/GalerkinToolkit.jl | 2 + src/assembly.jl | 132 ++++++++++++++------ src/interpolation.jl | 232 ++++++++++++++++-------------------- test/assembly_tests.jl | 8 +- test/interpolation_tests.jl | 18 +++ 6 files changed, 228 insertions(+), 166 deletions(-) diff --git a/Project.toml b/Project.toml index 2a7df090..ba9f1ef7 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,10 @@ authors = ["Francesc Verdugo and contributors"] version = "0.1.0" [deps] +BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" diff --git a/src/GalerkinToolkit.jl b/src/GalerkinToolkit.jl index b33e0301..eb6e0c01 100644 --- a/src/GalerkinToolkit.jl +++ b/src/GalerkinToolkit.jl @@ -9,6 +9,8 @@ import LinearAlgebra import ForwardDiff using Gmsh using PartitionedArrays +using BlockArrays +using FillArrays using Combinatorics using SparseArrays using Metis diff --git a/src/assembly.jl b/src/assembly.jl index 33197c58..4e3a978f 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -1,28 +1,90 @@ -# TODO monothinic matrix and vector for the moment -# add a strategy to build the matrix and vector +function monolithic_vector_assembly_strategy() + function init(dofs,::Type{T}) where T + n_total_dofs = length(dofs) + offsets = [0] + (;offsets,n_total_dofs,T) + end + function init(block_dofs::BRange,::Type{T}) where T + n_total_dofs = length(block_dofs) + offsets = blocklasts(block_dofs) .- map(length,blocks(block_dofs)) + (;offsets,n_total_dofs,T) + end + function scalar_type(setup) + setup.T + end + function counter(setup) + 0 + end + function count(n,setup,i,field_i) + n+1 + end + function allocate(n,setup) + Ti = Int32 + T = setup.T + I = zeros(Ti,n) + V = zeros(T,n) + (;I,V) + end + function set!(alloc,n,setup,v,i,field) + n += 1 + alloc.I[n] = i+setup.offsets[field] + alloc.V[n] = v + n + end + function compress(alloc,setup) + I = alloc.I + V = alloc.V + n_total_dofs = setup.n_total_dofs + vec = PartitionedArrays.dense_vector(I,V,n_total_dofs) + cache = nothing + (vec, cache) + end + function compress!(vec,cache,alloc,setup) + I = alloc.I + V = alloc.V + PartitionedArrays.dense_vector!(vec,I,V) + vec + end + (;init,scalar_type,counter,count,allocate,set!,compress,compress!) +end -function assemble_vector(f,space;kwargs...) +function assemble_vector(f,space,::Type{T};kwargs...) where T dim = 1 dv = GT.shape_functions(space,dim) integral = f(dv) - assemble_vector(integral,space;kwargs...) + assemble_vector(integral,space,T;kwargs...) end -function assemble_vector(integral::Number,space;reuse=false,Ti=Int32,T=Float64) +function assemble_vector(integral::Number,space,::Type{T}; + reuse = Val(false), + vector_strategy = monolithic_vector_assembly_strategy(), + ) where T @assert integral == 0 free_dofs = GT.free_dofs(space) - n = length(free_dofs) - zeros(T,n) + setup = vector_strategy.init(free_dofs,T) + counter = vector_strategy.counter(setup) + alloc = vector_strategy.allocate(counter,setup) + vec,cache = vector_strategy.compress(alloc,counter,setup) + if val_parameter(reuse) == false + vec + else + vec, (cache,alloc,setup,vector_strategy) + end end -function assemble_vector(integral::Integral,space;reuse=false,Ti=Int32,T=Float64) - state0 = (;integral,space,Ti,T) +function assemble_vector(integral::Integral,space,::Type{T}; + reuse = Val(false), + vector_strategy = monolithic_vector_assembly_strategy(), + ) where T + + setup = vector_strategy.init(GT.free_dofs(space),T) + state0 = (;integral,space,vector_strategy,setup) state1 = assemble_vector_count(state0) state2 = assemble_vector_allocate(state1) state3 = assemble_vector_fill(state2) result, cache = assemble_vector_compress(state3) - if reuse == false + if val_parameter(reuse) == false result else result, cache @@ -30,7 +92,7 @@ function assemble_vector(integral::Integral,space;reuse=false,Ti=Int32,T=Float64 end function assemble_vector_count(state) - (;integral,space) = state + (;integral,space,vector_strategy,setup) = state contributions = GT.contributions(integral) num_fields = GT.num_fields(space) fields = ntuple(identity,num_fields) @@ -41,7 +103,7 @@ function assemble_vector_count(state) GT.domain_glue(domain,GT.domain(space,field);strict=false) end end - n = 0 + counter = vector_strategy.counter(setup) function loop(glue::Nothing,field) end function loop(glue,field) @@ -51,8 +113,12 @@ function assemble_vector_count(state) for sface in 1:nsfaces faces = sface_to_faces[sface] for face in faces - ndofs = length(face_to_dofs[face]) - n += ndofs + dofs = face_to_dofs[face] + ndofs = length(dofs) + for idof in 1:ndofs + dof = dofs[idof] + counter = vector_strategy.count(counter,setup,dof,field) + end end end end @@ -60,21 +126,20 @@ function assemble_vector_count(state) domain, _ = domain_and_contribution map(loop,field_to_glue,fields) end - (;n,glues,fields,state...) + (;counter,glues,fields,state...) end function assemble_vector_allocate(state) - (;n,Ti,T) = state - I = zeros(Ti,n) - V = zeros(T,n) - (;I,V,state...) + (;counter,vector_strategy,setup) = state + alloc = vector_strategy.allocate(counter,setup) + (;alloc,state...) end function assemble_vector_fill(state) - (;integral,space,glues,I,V,fields) = state + (;integral,space,glues,fields,vector_strategy,alloc,setup) = state contributions = GT.contributions(integral) num_fields = GT.num_fields(space) - n = 0 + counter = vector_strategy.counter(setup) function loop(glue::Nothing,field,measure,contribution) end function loop(glue,field,measure,contribution) @@ -96,13 +161,13 @@ function assemble_vector_fill(state) s_qty = GT.topological_sort(expr_qty,(sface,face_around,idof,point)) s_npoints = GT.topological_sort(expr_npoints,(sface,)) expr = quote - function loop!(n,args,storage) - (;sface_to_faces,sface_to_faces_around,face_to_dofs,I,V) = args + function loop!(counter,args,storage) + (;sface_to_faces,sface_to_faces_around,face_to_dofs,vector_strategy,alloc,setup) = args $(unpack_storage(index.dict,:storage)) $(s_qty[1]) $(s_npoints[1]) nsfaces = length(sface_to_faces) - z = zero(eltype(V)) + z = zero(vector_strategy.scalar_type(setup)) for $sface in 1:nsfaces $(s_qty[2]) npoints = $(s_npoints[2]) @@ -118,20 +183,18 @@ function assemble_vector_fill(state) for $point in 1:npoints v += $(s_qty[5]) end - n += 1 dof = dofs[$idof] - I[n] = dof - V[n] = v + counter = vector_strategy.set!(alloc,counter,setup,v,dof,$field) end end end - n + counter end end loop! = eval(expr) - args = (;sface_to_faces,sface_to_faces_around,face_to_dofs,I,V) + args = (;sface_to_faces,sface_to_faces_around,face_to_dofs,vector_strategy,alloc,setup) storage = GT.storage(index) - n = invokelatest(loop!,n,args,storage) + counter = invokelatest(loop!,counter,args,storage) nothing end for (measure_and_contribution,field_to_glue) in zip(contributions,glues) @@ -144,12 +207,9 @@ function assemble_vector_fill(state) end function assemble_vector_compress(state) - (;I,V,space) = state - free_dofs = GT.free_dofs(space) - n = length(free_dofs) - vec = PartitionedArrays.dense_vector(I,V,n) - cache = (;I,V) - vec, cache + (;alloc,setup,vector_strategy) = state + vec, cache = vector_strategy.compress(alloc,setup) + vec, (cache,alloc,setup,vector_strategy) end function assemble_matrix(f,trial_space,test_space;kwargs...) diff --git a/src/interpolation.jl b/src/interpolation.jl index e6a3e06f..e8a121f9 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -139,12 +139,24 @@ Base.length(m::AbstractSpace) = num_fields(m) mesh(a::AbstractSpace) = GT.mesh(GT.domain(a)) +function free_dofs(a::AbstractSpace,field) + @assert field == 1 + free_dofs(a) +end + function free_dofs(a::AbstractSpace) - GT.dofs(free_values_strategy(a)) + nfree = length(first(free_and_dirichlet_dofs(a))) + Base.OneTo(nfree) +end + +function dirichlet_dofs(a::AbstractSpace,field) + @assert field == 1 + dirichlet_dofs(a) end function dirichlet_dofs(a::AbstractSpace) - GT.dofs(dirichlet_values_strategy(a)) + ndiri = length(last(free_and_dirichlet_dofs(a))) + Base.OneTo(ndiri) end # Single field by default @@ -695,47 +707,12 @@ function dual_basis(a::AbstractSpace,dim) u -> qty(compose(u,ϕ)) end -struct VectorStrategy{A,B,C,D} - dofs::A - allocate_values::B - restrict_values::C - prolongate_dof::D -end - -dofs(a::VectorStrategy) = a.dofs -allocate_values(a::VectorStrategy,args...) = a.allocate_values(args...) -restrict_values(a::VectorStrategy,args...) = a.restrict_values(args...) -prolongate_dof(a::VectorStrategy,args...) = a.prolongate_dof(args...) - -function monolithic_field_major_strategy(field_to_dofs) - offset = 0 - field_to_offset = map(field_to_dofs) do dofs - myoffset = offset - offset += length(dofs) - return myoffset - end - ndofs = offset - dofs = Base.OneTo(ndofs) - allocate_values(::Type{T}) where T = Vector{T}(undef,ndofs) - restrict_values(v,field) = view(v,field_to_dofs[field]) - prolongate_dof(dof,field) = dof + field_to_offset[field] - VectorStrategy(dofs,allocate_values,restrict_values,prolongate_dof) +function cartesian_product(spaces::AbstractSpace...) + CartesianProductSpace(spaces) end -function cartesian_product(spaces::AbstractSpace...; - free_values_strategy = GT.monolithic_field_major_strategy, - dirichlet_values_strategy = GT.monolithic_field_major_strategy, - ) - CartesianProductSpace( - spaces, - free_values_strategy(map(GT.free_dofs,spaces)), - dirichlet_values_strategy(map(GT.dirichlet_dofs,spaces)),) -end - -struct CartesianProductSpace{A,B,C} <: GT.AbstractSpace +struct CartesianProductSpace{A} <: GT.AbstractSpace spaces::A - free_values_strategy::B - dirichlet_values_strategy::C end function mesh(space::CartesianProductSpace) @@ -754,14 +731,6 @@ function LinearAlgebra.:×(a::AbstractSpace,b::CartesianProductSpace) cartesian_product(a,b.spaces...) end -function free_values_strategy(a::CartesianProductSpace) - a.free_values_strategy -end - -function dirichlet_values_strategy(a::CartesianProductSpace) - a.dirichlet_values_strategy -end - function components(a::CartesianProductSpace) a.spaces end @@ -792,45 +761,6 @@ function domain(a::CartesianProductSpace,field) GT.domain(GT.component(a,field)) end -#function num_face_dofs(a::CartesianProductSpace,dim) -# error("Casenot implemented (perhaps not needed in practice)") -#end -# -#function num_face_dofs(a::CartesianProductSpace,dim,field) -# f = component(a,field) -# num_face_dofs(f,dim,field) -# #terms = map(s->GT.num_face_dofs(s,dim),a.spaces) -# #index -> begin -# # field = index.field_per_dim[dim] -# # @assert field !== nothing -# # index2 = replace_field_per_dim(index,nothing) -# # terms[field](index2) -# #end -#end - -## TODO face_dofs vs dof_map very confusing. Remove dof_map? -#function dof_map(a::CartesianProductSpace,dim) -# error("Casenot implemented (perhaps not needed in practice)") -#end -# -#function dof_map(a::CartesianProductSpace,dim,field) -# strategy = GT.free_values_strategy(a) -# f = component(a,field) -# term = dof_map(f,dim,field) -# index -> begin -# dof_term = term(index) -# dict = index.dict -# field_sym = get!(dict,field,gensym("field")) -# strategy_sym = get!(dict,strategy,gensym("strategy")) -# @term begin -# dof = $dof_term -# prolongate_dof($strategy_sym,dof,$field_sym) -# end -# -# end -#end -# - function face_dofs(a::CartesianProductSpace) error("Not implemented, not needed in practice") end @@ -839,6 +769,30 @@ function face_dofs(a::CartesianProductSpace,field) face_dofs(a.spaces[field]) end +function free_dofs(a::CartesianProductSpace) + nfields = GT.num_fields(a) + map(1:nfields) do field + free_dofs(a,field) + end |> BRange +end + +function free_dofs(a::CartesianProductSpace,field) + f = component(a,field) + free_dofs(f) +end + +function dirichlet_dofs(a::CartesianProductSpace) + nfields = GT.num_fields(a) + map(1:nfields) do field + dirichlet_dofs(a,field) + end |> BRange +end + +function dirichlet_dofs(a::CartesianProductSpace,field) + f = component(a,field) + dirichlet_dofs(f) +end + function shape_functions(a::CartesianProductSpace,dim) fields = ntuple(identity,GT.num_fields(a)) map(fields) do field @@ -892,9 +846,31 @@ Base.iterate(m::DiscreteField,state) = iterate(components(m),state) Base.getindex(m::DiscreteField,field::Integer) = component(m,field) Base.length(m::DiscreteField) = num_fields(m) +function allocate_values(::Type{T},dofs) where T + n = length(dofs) + Vector{T}(undef,n) +end + +function allocate_values(::Type{T},dofs::BRange) where T + map(dofs.blocks) do dofs_i + allocate_values(T,dofs_i) + end |> BVector +end + +function constant_values(f,dofs) + n = length(dofs) + Fill(f,n) +end + +function constant_values(f,dofs::BRange) + map(dofs.blocks) do dofs_i + constant_values(f,dofs_i) + end |> BVector +end + function undef_field(::Type{T},space::AbstractSpace) where T - free_values = allocate_values(free_values_strategy(space),T) - dirichlet_values = allocate_values(dirichlet_values_strategy(space),T) + free_values = allocate_values(T,GT.free_dofs(space)) + dirichlet_values = allocate_values(T,GT.dirichlet_dofs(space)) discrete_field(space,free_values,dirichlet_values) end @@ -915,6 +891,21 @@ end # TODO rand_field +function undef_dirichlet_field(::Type{T},space::AbstractSpace) where T + free_values = constant_values(zero(T),GT.free_dofs(space)) + dirichlet_values = allocate_values(T,GT.dirichlet_dofs(space)) + discrete_field(space,free_values,dirichlet_values) +end + +function zero_dirichlet_field(::Type{T},space::AbstractSpace) where T + uh = undef_dirichlet_field(T,space) + fill!(dirichlet_values(uh),zero(T)) +end + +function dirichlet_field(::Type{T},space::AbstractSpace) where T + zero_dirichlet_field(T,space) +end + function num_fields(u::DiscreteField) num_fields(GT.space(u)) end @@ -928,8 +919,8 @@ end function component(u::DiscreteField,field::Integer) space = GT.space(u) space_field = component(space,field) - free_values_field = restrict_values(free_values_strategy(space),free_values(u),field) - dirichlet_values_field = restrict_values(dirichlet_values_strategy(space),dirichlet_values(u),field) + free_values_field = view(free_values(u),Block(field)) + dirichlet_values_field = view(dirichlet_values(u),Block(field)) discrete_field(space_field,free_values_field,dirichlet_values_field) end @@ -1013,19 +1004,38 @@ function interpolate!(f,u::DiscreteField) interpolate!(f,u,nothing) end +function interpolate!(f,u::DiscreteField,field) + interpolate!(f,u,nothing,field) +end + function interpolate_free!(f,u::DiscreteField) interpolate!(f,u,FREE) end +function interpolate_free!(f,u::DiscreteField,field) + interpolate!(f,u,FREE,field) +end + function interpolate_dirichlet!(f,u::DiscreteField) # TODO for dirichlet we perhaps want to allow integrating on boundaries? interpolate!(f,u,DIRICHLET) end +function interpolate_dirichlet!(f,u::DiscreteField,field) + # TODO for dirichlet we perhaps want to allow integrating on boundaries? + interpolate!(f,u,DIRICHLET,field) +end + function interpolate!(f,u::DiscreteField,free_or_diri::Union{Nothing,FreeOrDirichlet}) interpolate_impl!(f,u,free_or_diri) end +function interpolate!(f,u::DiscreteField,free_or_diri::Union{Nothing,FreeOrDirichlet},field) + ui = component(u,field) + interpolate_impl!(f,ui,free_or_diri) + u +end + function interpolate_impl!(f,u,free_or_diri;location=1) free_vals = GT.free_values(u) diri_vals = GT.dirichlet_values(u) @@ -1084,20 +1094,14 @@ end # This space is not suitable for assembling problems in general, as there might be gaps in the dofs -# TODO think about free_values_strategy and dirichlet_values_strategy function iso_parametric_space(dom::ReferenceDomain; - dirichlet_boundary=nothing, - free_values_strategy = GT.monolithic_field_major_strategy, - dirichlet_values_strategy = GT.monolithic_field_major_strategy, - ) - IsoParametricSpace(dom,dirichlet_boundary,free_values_strategy,dirichlet_values_strategy) + dirichlet_boundary=nothing) + IsoParametricSpace(dom,dirichlet_boundary) end -struct IsoParametricSpace{A,B,C,D} <: AbstractSpace +struct IsoParametricSpace{A,B} <: AbstractSpace domain::A dirichlet_boundary::B - free_values_strategy::C - dirichlet_values_strategy::D end function free_and_dirichlet_nodes(V::IsoParametricSpace) @@ -1172,16 +1176,6 @@ function face_reference_id(V::IsoParametricSpace) view(face_refid,faces) end -function free_values_strategy(a::IsoParametricSpace) - n = length(first(free_and_dirichlet_nodes(a))) - a.free_values_strategy( (Base.OneTo(n),) ) -end - -function dirichlet_values_strategy(a::IsoParametricSpace) - n = length(last(free_and_dirichlet_nodes(a))) - a.dirichlet_values_strategy( (Base.OneTo(n),) ) -end - function domain(a::IsoParametricSpace) a.domain end @@ -1190,8 +1184,6 @@ end function lagrange_space(domain,order; conformity = :default, dirichlet_boundary=nothing, - free_values_strategy = GT.monolithic_field_major_strategy, - dirichlet_values_strategy = GT.monolithic_field_major_strategy, space=nothing, major=:component, shape=SCALAR_SHAPE) @@ -1203,21 +1195,17 @@ function lagrange_space(domain,order; order, conformity, dirichlet_boundary, - free_values_strategy, - dirichlet_values_strategy, space, major, shape) end -struct LagrangeSpace{A,B,C,D,E,F,G,H} <: AbstractSpace +struct LagrangeSpace{A,B,C,F,G,H} <: AbstractSpace domain::A order::B conformity::Symbol dirichlet_boundary::C - free_values_strategy::D - dirichlet_values_strategy::E space::F # TODO rename this one? major::G shape::H @@ -1260,13 +1248,3 @@ function reference_fes(space::LagrangeSpace) ctype_to_reffe end -function free_values_strategy(a::LagrangeSpace) - n = length(first(free_and_dirichlet_dofs(a))) - a.free_values_strategy( (Base.OneTo(n),) ) -end - -function dirichlet_values_strategy(a::LagrangeSpace) - n = length(last(free_and_dirichlet_dofs(a))) - a.dirichlet_values_strategy( (Base.OneTo(n),) ) -end - diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index e9b56e21..41bdc913 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -73,7 +73,7 @@ function l(v) end end -b = GT.assemble_vector(l,V) +b = GT.assemble_vector(l,V,Float64) function l(v) ∫(dΛref) do p @@ -82,7 +82,7 @@ function l(v) end end -b = GT.assemble_vector(l,V) +b = GT.assemble_vector(l,V,Float64) @test sum(b)+1 ≈ 1 V² = V × V @@ -98,7 +98,9 @@ function l((v1,v2)) end end -b = GT.assemble_vector(l,V²) +b = GT.assemble_vector(l,V²,Float64) + +xxxx function a(u,v) ∫(dΩref) do q diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl index 37f3f40e..47e2bf87 100644 --- a/test/interpolation_tests.jl +++ b/test/interpolation_tests.jl @@ -5,6 +5,7 @@ using GalerkinToolkit: ∫, × using Test import ForwardDiff using LinearAlgebra +using BlockArrays D = 2 order = 3 @@ -65,8 +66,22 @@ GT.interpolate_dirichlet!(uref,v2) Y = V×V y = GT.zero_field(Float64,Y) +display(GT.free_values(y)) +display(GT.free_dofs(Y)) y1, y2 = y +@test GT.free_values(y1) === blocks(GT.free_values(y))[1] +@test GT.free_values(y2) === blocks(GT.free_values(y))[2] + +GT.interpolate!(uref,y1) +GT.interpolate!(uref,y,1) + +GT.interpolate_free!(uref,y1) +GT.interpolate_free!(uref,y,1) + +GT.interpolate_dirichlet!(uref,y1) +GT.interpolate_dirichlet!(uref,y,1) + order = 3 V = GT.lagrange_space(Ωref,order) GT.face_dofs(V) @@ -98,6 +113,9 @@ end el2 = ∫( q->abs2(eh(q))*dV(q), dΩref) |> sum |> sqrt @test el2 < tol +uhd = GT.dirichlet_field(Float64,V) +GT.interpolate_dirichlet!(uref,uh) + # TODO #udiri = GT.analytical_field(sum,Γdiri) #uh = GT.zero_field(Float64,V) From 3cbb84a2808174e8aa90f21bbe1a50cdd2fb8140 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 28 Aug 2024 18:11:11 +0200 Subject: [PATCH 15/20] Implemented matrix assembly strategy --- src/assembly.jl | 143 +++++++++++++++++++++++++++++------------ test/assembly_tests.jl | 7 +- 2 files changed, 105 insertions(+), 45 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 4e3a978f..4dfb0740 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -49,6 +49,64 @@ function monolithic_vector_assembly_strategy() (;init,scalar_type,counter,count,allocate,set!,compress,compress!) end +function monolithic_matrix_assembly_strategy(;matrix_type=nothing) + function init(dofs_test,dofs_trial,::Type{T}) where T + n_total_rows = length(dofs_test) + n_total_cols = length(dofs_trial) + offsets_rows = [0] + offsets_cols = [0] + (;offsets_rows,offsets_cols,n_total_rows,n_total_cols,T) + end + function init(dofs_test::BRange,dofs_trial,::Type{T}) where T + n_total_rows = length(dofs_test) + n_total_cols = length(dofs_trial) + offsets_rows = blocklasts(dofs_test) .- map(length,blocks(dofs_test)) + offsets_cols = blocklasts(dofs_trial) .- map(length,blocks(dofs_trial)) + (;offsets_rows,offsets_cols,n_total_rows,n_total_cols,T) + end + function scalar_type(setup) + setup.T + end + function counter(setup) + 0 + end + function count(n,setup,i,j,field_i,field_j) + n+1 + end + function allocate(n,setup) + Ti = Int32 + T = setup.T + I = zeros(Ti,n) + J = zeros(Ti,n) + V = zeros(T,n) + (;I,J,V) + end + function set!(alloc,n,setup,v,i,j,field_i,field_j) + n += 1 + alloc.I[n] = i+setup.offsets_rows[field_i] + alloc.J[n] = j+setup.offsets_cols[field_j] + alloc.V[n] = v + n + end + function compress(alloc,setup) + I = alloc.I + J = alloc.J + V = alloc.V + n_total_rows = setup.n_total_rows + n_total_cols = setup.n_total_cols + T = setup.T + M = matrix_type === nothing ? SparseMatrixCSC{T,Int} : matrix_type + A,cache = sparse_matrix(M,I,J,V,n_total_rows,n_total_cols;reuse=Val(true)) + (A, cache) + end + function compress!(A,cache,alloc,setup) + V = alloc.V + sparse_matrix!(A,V,cache) + A + end + (;init,scalar_type,counter,count,allocate,set!,compress,compress!) +end + function assemble_vector(f,space,::Type{T};kwargs...) where T dim = 1 dv = GT.shape_functions(space,dim) @@ -65,7 +123,7 @@ function assemble_vector(integral::Number,space,::Type{T}; setup = vector_strategy.init(free_dofs,T) counter = vector_strategy.counter(setup) alloc = vector_strategy.allocate(counter,setup) - vec,cache = vector_strategy.compress(alloc,counter,setup) + vec,cache = vector_strategy.compress(alloc,setup) if val_parameter(reuse) == false vec else @@ -162,7 +220,7 @@ function assemble_vector_fill(state) s_npoints = GT.topological_sort(expr_npoints,(sface,)) expr = quote function loop!(counter,args,storage) - (;sface_to_faces,sface_to_faces_around,face_to_dofs,vector_strategy,alloc,setup) = args + (;sface_to_faces,sface_to_faces_around,face_to_dofs,field,vector_strategy,alloc,setup) = args $(unpack_storage(index.dict,:storage)) $(s_qty[1]) $(s_npoints[1]) @@ -184,7 +242,7 @@ function assemble_vector_fill(state) v += $(s_qty[5]) end dof = dofs[$idof] - counter = vector_strategy.set!(alloc,counter,setup,v,dof,$field) + counter = vector_strategy.set!(alloc,counter,setup,v,dof,field) end end end @@ -192,7 +250,7 @@ function assemble_vector_fill(state) end end loop! = eval(expr) - args = (;sface_to_faces,sface_to_faces_around,face_to_dofs,vector_strategy,alloc,setup) + args = (;sface_to_faces,sface_to_faces_around,face_to_dofs,field,vector_strategy,alloc,setup) storage = GT.storage(index) counter = invokelatest(loop!,counter,args,storage) nothing @@ -212,17 +270,21 @@ function assemble_vector_compress(state) vec, (cache,alloc,setup,vector_strategy) end -function assemble_matrix(f,trial_space,test_space;kwargs...) +function assemble_matrix(f,trial_space,test_space,::Type{T};kwargs...) where T test_dim = 1 trial_dim = 2 dv = GT.shape_functions(test_space,test_dim) du = GT.shape_functions(trial_space,trial_dim) integral = f(du,dv) - assemble_matrix(integral,trial_space,test_space;kwargs...) + assemble_matrix(integral,trial_space,test_space,T;kwargs...) end -function assemble_matrix(integral::Integral,trial_space,test_space;reuse=false,Ti=Int32,T=Float64) - state0 = (;integral,test_space,trial_space,Ti,T) +function assemble_matrix(integral::Integral,trial_space,test_space,::Type{T}; + reuse=false, + matrix_strategy = monolithic_matrix_assembly_strategy(), + ) where T + setup = matrix_strategy.init(free_dofs(test_space),free_dofs(trial_space),T) + state0 = (;integral,test_space,trial_space,matrix_strategy,setup) state1 = assemble_matrix_count(state0) state2 = assemble_matrix_allocate(state1) state3 = assemble_matrix_fill(state2) @@ -235,7 +297,7 @@ function assemble_matrix(integral::Integral,trial_space,test_space;reuse=false,T end function assemble_matrix_count(state) - (;integral,test_space,trial_space) = state + (;integral,test_space,trial_space,matrix_strategy,setup) = state contributions = GT.contributions(integral) num_fields_test = GT.num_fields(test_space) num_fields_trial = GT.num_fields(trial_space) @@ -255,13 +317,14 @@ function assemble_matrix_count(state) GT.domain_glue(domain,GT.domain(trial_space,field),strict=false) end end - n = 0 + counter = matrix_strategy.counter(setup) function loop(glue_test,glue_trial,field_per_dim) sface_to_faces_test, = target_face(glue_test) sface_to_faces_trial, = target_face(glue_trial) face_to_dofs_test = face_dofs(test_space,field_per_dim[1]) face_to_dofs_trial = face_dofs(trial_space,field_per_dim[2]) nsfaces = length(sface_to_faces_test) + field_test,field_trial = field_per_dim for sface in 1:nsfaces faces_test = sface_to_faces_test[sface] faces_trial = sface_to_faces_trial[sface] @@ -271,7 +334,13 @@ function assemble_matrix_count(state) for face_trial in faces_trial dofs_trial = face_to_dofs_trial[face_trial] ndofs_trial = length(dofs_trial) - n += ndofs_test*ndofs_trial + for idof_test in 1:ndofs_test + dof_test = dofs_test[idof_test] + for idof_trial in 1:ndofs_trial + dof_trial = dofs_trial[idof_trial] + counter = matrix_strategy.count(counter,setup,dof_test,dof_trial,field_test,field_trial) + end + end end end end @@ -290,28 +359,27 @@ function assemble_matrix_count(state) end end end - (;n,glues_test,glues_trial,fields_test,fields_trial,state...) + (;counter,glues_test,glues_trial,fields_test,fields_trial,state...) end function assemble_matrix_allocate(state) - (;n,Ti,T) = state - I = zeros(Ti,n) - J = zeros(Ti,n) - V = zeros(T,n) - (;I,J,V,state...) + (;matrix_strategy,counter,setup) = state + alloc = matrix_strategy.allocate(counter,setup) + (;alloc,state...) end function assemble_matrix_fill(state) - (;integral,test_space,trial_space,fields_test,fields_trial,I,J,V,glues_test,glues_trial) = state + (;integral,test_space,trial_space,fields_test,fields_trial,alloc,glues_test,glues_trial,matrix_strategy,setup) = state contributions = GT.contributions(integral) num_fields_test = GT.num_fields(test_space) num_fields_trial = GT.num_fields(trial_space) - n = 0 + counter = matrix_strategy.counter(setup) function loop(glue_test,glue_trial,field_per_dim,measure,contribution) sface_to_faces_test, _,sface_to_faces_around_test = target_face(glue_test) sface_to_faces_trial, _,sface_to_faces_around_trial = target_face(glue_trial) face_to_dofs_test = face_dofs(test_space,field_per_dim[1]) face_to_dofs_trial = face_dofs(trial_space,field_per_dim[2]) + field_test,field_trial = field_per_dim term_qty = GT.term(contribution) term_npoints = GT.num_points(measure) sface = :sface @@ -329,15 +397,16 @@ function assemble_matrix_fill(state) s_qty = GT.topological_sort(expr_qty,(sface,face_around_test,face_around_trial,idof_test,idof_trial,point)) s_npoints = GT.topological_sort(expr_npoints,(sface,)) expr = quote - function loop!(n,args,storage) + function loop!(counter,args,storage) (;sface_to_faces_test,sface_to_faces_around_test, sface_to_faces_trial,sface_to_faces_around_trial, - face_to_dofs_test,face_to_dofs_trial,I,J,V) = args + face_to_dofs_test,face_to_dofs_trial,field_test,field_trial, + alloc,matrix_strategy,setup) = args $(unpack_storage(index.dict,:storage)) $(s_qty[1]) $(s_npoints[1]) nsfaces = length(sface_to_faces_test) - z = zero(eltype(V)) + z = zero(matrix_strategy.scalar_type(setup)) for $sface in 1:nsfaces $(s_qty[2]) npoints = $(s_npoints[2]) @@ -363,24 +432,22 @@ function assemble_matrix_fill(state) for $point in 1:npoints v += $(s_qty[7]) end - n += 1 - I[n] = dof_test - J[n] = dof_trial - V[n] = v + counter = matrix_strategy.set!(alloc,counter,setup,v,dof_test,dof_trial,field_test,field_trial) end end end end end - n + counter end end loop! = eval(expr) storage = GT.storage(index) args = (;sface_to_faces_test,sface_to_faces_around_test, sface_to_faces_trial,sface_to_faces_around_trial, - face_to_dofs_test,face_to_dofs_trial,I,J,V) - n = invokelatest(loop!,n,args,storage) + face_to_dofs_test,face_to_dofs_trial,field_test,field_trial, + alloc,matrix_strategy,setup) + counter = invokelatest(loop!,counter,args,storage) nothing end for (measure_and_contribution,field_to_glue_test,field_to_glue_trial) in zip(contributions,glues_test,glues_trial) @@ -401,25 +468,19 @@ function assemble_matrix_fill(state) end function assemble_matrix_compress(state) - (;I,J,V,test_space,trial_space) = state - free_dofs_test = GT.free_dofs(test_space) - free_dofs_trial = GT.free_dofs(trial_space) - m = length(free_dofs_test) - n = length(free_dofs_trial) - vec = PartitionedArrays.sparse_matrix(I,J,V,m,n) - cache = (;I,J,V) - vec, cache + (;alloc,matrix_strategy,setup) = state + A, cache = matrix_strategy.compress(alloc,setup) + A, (cache,alloc,setup,matrix_strategy) end - function linear_problem(uh,a,l,V=GT.space(uh)) U = GT.space(uh) x = free_values(uh) fill!(x,0) T = eltype(x) - A = assemble_matrix(a,U,V;T) - b = assemble_vector(l,V;T) - d = assemble_vector(v->a(uh,v),V) + A = assemble_matrix(a,U,V,T) + b = assemble_vector(l,V,T) + d = assemble_vector(v->a(uh,v),V,T) b .= b .- d x,A,b end diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index 41bdc913..a2c1e718 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -100,8 +100,6 @@ end b = GT.assemble_vector(l,V²,Float64) -xxxx - function a(u,v) ∫(dΩref) do q J = ForwardDiff.jacobian(ϕ,q) @@ -118,7 +116,7 @@ function a(u,v) end end -A = GT.assemble_matrix(a,V,V) +A = GT.assemble_matrix(a,V,V,Float64) function a((u1,u2),(v1,v2)) ∫(dΩref) do q @@ -131,7 +129,7 @@ function a((u1,u2),(v1,v2)) end end -A = GT.assemble_matrix(a,V²,V²) +A = GT.assemble_matrix(a,V²,V²,Float64) function dV(ϕ,q) J = ForwardDiff.jacobian(ϕ,q) @@ -160,6 +158,7 @@ a(u,v) = ∫( q->∇(u,ϕ,q)⋅∇(v,ϕ,q)*dV(ϕ,q), dΩref) l(v) = 0 x,A,b = GT.linear_problem(uh,a,l) +display(A) x .= A\b # Poisson solve (reference domain) From 3f98c52d58e95faf64954773ffe7a20637eaef6c Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 28 Aug 2024 18:27:32 +0200 Subject: [PATCH 16/20] Introduced solution_field --- src/assembly.jl | 29 +++++++++++++++++++++++------ src/interpolation.jl | 1 + test/assembly_tests.jl | 18 ++++++++++-------- 3 files changed, 34 insertions(+), 14 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 4dfb0740..0af88832 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -473,15 +473,32 @@ function assemble_matrix_compress(state) A, (cache,alloc,setup,matrix_strategy) end -function linear_problem(uh,a,l,V=GT.space(uh)) - U = GT.space(uh) - x = free_values(uh) - fill!(x,0) - T = eltype(x) +function linear_problem(uhd::DiscreteField,a,l,V=GT.space(uhd)) + U = GT.space(uhd) + T = eltype(dirichlet_values(uhd)) A = assemble_matrix(a,U,V,T) b = assemble_vector(l,V,T) - d = assemble_vector(v->a(uh,v),V,T) + d = assemble_vector(v->a(uhd,v),V,T) b .= b .- d + x = similar(b,axes(A,2)) + x,A,b +end + +function linear_problem(::Type{T},U::AbstractSpace,a,l,V=U) where T + A = assemble_matrix(a,U,V,T) + b = assemble_vector(l,V,T) + x = similar(b,axes(A,2)) x,A,b end +function solution_field(uhd::DiscreteField,x) + diri_values = dirichlet_values(uhd) + U = GT.space(uhd) + discrete_field(U,x,diri_values) +end + +function solution_field(U::AbstractSpace,x) + T = eltype(x) + uhd = zero_dirichlet_field(T,U) + solution_field(uhd,x) +end diff --git a/src/interpolation.jl b/src/interpolation.jl index e8a121f9..33a94a63 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -900,6 +900,7 @@ end function zero_dirichlet_field(::Type{T},space::AbstractSpace) where T uh = undef_dirichlet_field(T,space) fill!(dirichlet_values(uh),zero(T)) + uh end function dirichlet_field(::Type{T},space::AbstractSpace) where T diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index a2c1e718..2822852c 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -143,10 +143,10 @@ f = GT.analytical_field(sum,Ω) l(v) = ∫( q->f(ϕ(q))*v(q)*dV(ϕ,q), dΩref) V = GT.iso_parametric_space(Ωref) -uh = GT.zero_field(Float64,V) -x,A,b = GT.linear_problem(uh,a,l) +x,A,b = GT.linear_problem(Float64,V,a,l) x .= A\b +uh = GT.solution_field(V,x) function ∇(u,phi,q) J = ForwardDiff.jacobian(phi,q) @@ -181,10 +181,10 @@ order = 3 V = GT.lagrange_space(Ωref,order;dirichlet_boundary=Γdiri) u = GT.analytical_field(sum,Ω) -uh = GT.zero_field(Float64,V) +uhd = GT.dirichlet_field(Float64,V) # TODO #GT.interpolate_dirichlet!(q->u(ϕ(q)),uh) -GT.interpolate_dirichlet!(u∘ϕ,uh) +GT.interpolate_dirichlet!(u∘ϕ,uhd) function ∇(u,q) J = ForwardDiff.jacobian(ϕ,q) @@ -203,8 +203,9 @@ dΩref = GT.measure(Ωref,degree) a(u,v) = ∫( q->∇(u,q)⋅∇(v,q)*dV(q), dΩref) l(v) = 0 -x,A,b = GT.linear_problem(uh,a,l) +x,A,b = GT.linear_problem(uhd,a,l) x .= A\b +uh = GT.solution_field(uhd,x) # TODO # Functions like this ones should @@ -223,8 +224,8 @@ eh1 = ∫( q->∇eh(q)⋅∇eh(q)*dV(q), dΩref) |> sum |> sqrt V = GT.lagrange_space(Ω,order;dirichlet_boundary=Γdiri) -uh = GT.zero_field(Float64,V) -GT.interpolate_dirichlet!(u,uh) +uhd = GT.dirichlet_field(Float64,V) +GT.interpolate_dirichlet!(u,uhd) dΩ = GT.measure(Ω,degree) @@ -233,8 +234,9 @@ dΩ = GT.measure(Ω,degree) a(u,v) = ∫( q->∇(u,q)⋅∇(v,q), dΩ) l(v) = 0 -x,A,b = GT.linear_problem(uh,a,l) +x,A,b = GT.linear_problem(uhd,a,l) x .= A\b +uh = GT.solution_field(uhd,x) eh(q) = u(q) - uh(q) ∇eh(q) = ∇(u,q) - ∇(uh,q) From 7f0de2301e0ae380aa31173efdae462c1c161f7b Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Wed, 28 Aug 2024 20:03:48 +0200 Subject: [PATCH 17/20] Bugfix in solution_field --- src/assembly.jl | 32 ++++++++++++++++++++++++++------ test/assembly_tests.jl | 7 +++++++ 2 files changed, 33 insertions(+), 6 deletions(-) diff --git a/src/assembly.jl b/src/assembly.jl index 0af88832..ef8dfa6b 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -491,14 +491,34 @@ function linear_problem(::Type{T},U::AbstractSpace,a,l,V=U) where T x,A,b end -function solution_field(uhd::DiscreteField,x) - diri_values = dirichlet_values(uhd) - U = GT.space(uhd) - discrete_field(U,x,diri_values) -end - function solution_field(U::AbstractSpace,x) T = eltype(x) uhd = zero_dirichlet_field(T,U) solution_field(uhd,x) end + +function solution_field(uhd::DiscreteField,x) + diri_vals = dirichlet_values(uhd) + U = GT.space(uhd) + free_vals = free_values_from_solution(x,free_dofs(U)) + discrete_field(U,free_vals,diri_vals) +end + +function free_values_from_solution(x,dofs) + x +end + +function free_values_from_solution(x,dofs::BRange) + nfields = blocklength(dofs) + map(1:nfields) do field + pend = blocklasts(dofs)[field] + pini = 1 + pend - length(blocks(dofs)[field]) + view(x,pini:pend) + end |> BVector +end + +function free_values_from_solution(x::BVector,dofs::BRange) + x +end + + diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index 2822852c..b71f0bc1 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -131,6 +131,13 @@ end A = GT.assemble_matrix(a,V²,V²,Float64) +x = similar(b,axes(A,2)) +fill!(x,0) +uh = GT.solution_field(V²,x) +uh1,uh2 = uh +fill!(GT.free_values(uh2),1) +@test x[end] == 1 + function dV(ϕ,q) J = ForwardDiff.jacobian(ϕ,q) abs(det(J)) From 67b70133cba7a61904b937a4b46164125727fb71 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 29 Aug 2024 07:50:50 +0200 Subject: [PATCH 18/20] Most of Poisson example working (except for multipliers) --- GalerkinToolkitExamples/src/poisson.jl | 38 +++++++++++-------- GalerkinToolkitExamples/test/poisson_tests.jl | 2 +- Project.toml | 2 +- src/domain.jl | 2 +- 4 files changed, 26 insertions(+), 18 deletions(-) diff --git a/GalerkinToolkitExamples/src/poisson.jl b/GalerkinToolkitExamples/src/poisson.jl index 1bb26880..5e860cf2 100644 --- a/GalerkinToolkitExamples/src/poisson.jl +++ b/GalerkinToolkitExamples/src/poisson.jl @@ -34,8 +34,8 @@ laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x)) Δ(u,x) = Δ(u)(x) ∇(u,x) = ∇(u)(x) -mean(u,x) = 0.5*(u(x)[1]+u(x)[2]) -jump(u,n,x) = u(x)[1]*n[1](x) + u(x)[2]*n[2](x) +mean(u,x) = 0.5*(u(x)[+]+u(x)[-]) +jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x) function main_automatic(params) timer = params[:timer] @@ -74,21 +74,18 @@ function main_automatic(params) if params[:dirichlet_method] === :strong V = GT.lagrange_space(Ω,interpolation_degree;conformity,dirichlet_boundary=Γd) - uh = GT.zero_field(Float64,V) - GT.interpolate_dirichlet!(u,uh) + uhd = GT.dirichlet_field(Float64,V) + GT.interpolate_dirichlet!(u,uhd) else n_Γd = GT.unit_normal(Γd,Ω) h_Γd = GT.face_diameter_field(Γd) dΓd = GT.measure(Γd,integration_degree) V = GT.lagrange_space(Ω,interpolation_degree;conformity) - uh = GT.zero_field(Float64,V) end if params[:dirichlet_method] === :multipliers Q = GT.lagrange_space(Γd,interpolation_degree-1;conformity=:L2) VxQ = V × Q - uh_qh = GT.zero_field(Float64,VxQ) - uh, qh = uh_qh end function a(u,v) @@ -123,21 +120,32 @@ function main_automatic(params) r += ∫(x->u(x)*q(x), dΓd) r end - Uh = uh_qh - else - A = a - L = l - Uh = uh end @timeit timer "assembly" begin - # TODO give a hint when dirichlet BCS are homogeneous or not present - x,A,b = GT.linear_problem(Uh,A,L) + if params[:dirichlet_method] === :strong + # TODO give a hint when dirichlet BCS are homogeneous + x,K,b = GT.linear_problem(uhd,a,l) + elseif params[:dirichlet_method] === :multipliers + x,K,b = GT.linear_problem(Float64,VxQ,A,L) + else + x,K,b = GT.linear_problem(Float64,V,a,l) + end end @timeit timer "solver" begin - P = ps.setup(params[:solver],x,A,b) + P = ps.setup(params[:solver],x,K,b) + # TODO the solver should tell if the initial guess will be actually used + fill!(x,0) # set initial guess ps.solve!(x,P,b) + if params[:dirichlet_method] === :strong + uh = GT.solution_field(uhd,x) + elseif params[:dirichlet_method] === :multipliers + uh_ph = GT.solution_field(VxQ,x) + uh, ph = uh_ph + else + uh = GT.solution_field(V,x) + end end @timeit timer "error_norms" begin diff --git a/GalerkinToolkitExamples/test/poisson_tests.jl b/GalerkinToolkitExamples/test/poisson_tests.jl index 2fad66e5..2f337c77 100644 --- a/GalerkinToolkitExamples/test/poisson_tests.jl +++ b/GalerkinToolkitExamples/test/poisson_tests.jl @@ -11,7 +11,7 @@ tol = 1.0e-10 n = 2 mesh = GT.cartesian_mesh((0,2,0,2),(n,n)) for discretization_method in (:interior_penalty,:continuous_galerkin) - for dirichlet_method in (:multipliers,:nitsche,:strong) + for dirichlet_method in (:nitsche,:strong) #(:multipliers,:nitsche,:strong) params = Dict{Symbol,Any}() params[:mesh] = mesh params[:dirichlet_tags] = ["1-face-1","1-face-3","1-face-4"] diff --git a/Project.toml b/Project.toml index ba9f1ef7..5542a0dc 100644 --- a/Project.toml +++ b/Project.toml @@ -28,7 +28,7 @@ Gmsh = "0.3" IterativeSolvers = "0.9" MPI = "0.20" Metis = "1" -PartitionedArrays = "0.4, 0.5" +PartitionedArrays = "0.5.4" StaticArrays = "1" WriteVTK = "1" julia = "1" diff --git a/src/domain.jl b/src/domain.jl index b15a3028..c6149ef3 100644 --- a/src/domain.jl +++ b/src/domain.jl @@ -589,7 +589,7 @@ function analytical_field(f,dom::AbstractDomain) constant_quantity(f,dom) end -function face_constant_field_impl(date,face::Integer) +function face_constant_field_impl(data,face::Integer) x->data[face] end From f52262becb69f5380fc8c03753b58fd6d1c62b65 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 29 Aug 2024 16:09:12 +0200 Subject: [PATCH 19/20] Fixing examples --- GalerkinToolkitExamples/test/poisson_tests.jl | 2 +- src/assembly.jl | 24 ++++++++++++++++--- src/interpolation.jl | 14 +++++++++-- 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/GalerkinToolkitExamples/test/poisson_tests.jl b/GalerkinToolkitExamples/test/poisson_tests.jl index 2f337c77..2fad66e5 100644 --- a/GalerkinToolkitExamples/test/poisson_tests.jl +++ b/GalerkinToolkitExamples/test/poisson_tests.jl @@ -11,7 +11,7 @@ tol = 1.0e-10 n = 2 mesh = GT.cartesian_mesh((0,2,0,2),(n,n)) for discretization_method in (:interior_penalty,:continuous_galerkin) - for dirichlet_method in (:nitsche,:strong) #(:multipliers,:nitsche,:strong) + for dirichlet_method in (:multipliers,:nitsche,:strong) params = Dict{Symbol,Any}() params[:mesh] = mesh params[:dirichlet_tags] = ["1-face-1","1-face-3","1-face-4"] diff --git a/src/assembly.jl b/src/assembly.jl index ef8dfa6b..ab795d7a 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -171,6 +171,9 @@ function assemble_vector_count(state) for sface in 1:nsfaces faces = sface_to_faces[sface] for face in faces + if face == 0 + continue + end dofs = face_to_dofs[face] ndofs = length(dofs) for idof in 1:ndofs @@ -232,6 +235,9 @@ function assemble_vector_fill(state) faces = sface_to_faces[sface] faces_around = sface_to_faces_around[sface] for ($face_around,face) in zip(faces_around,faces) + if face == 0 + continue + end $(s_qty[3]) dofs = face_to_dofs[face] ndofs = length(dofs) @@ -329,9 +335,15 @@ function assemble_matrix_count(state) faces_test = sface_to_faces_test[sface] faces_trial = sface_to_faces_trial[sface] for face_test in faces_test + if face_test == 0 + continue + end dofs_test = face_to_dofs_test[face_test] ndofs_test = length(dofs_test) for face_trial in faces_trial + if face_trial == 0 + continue + end dofs_trial = face_to_dofs_trial[face_trial] ndofs_trial = length(dofs_trial) for idof_test in 1:ndofs_test @@ -375,11 +387,11 @@ function assemble_matrix_fill(state) num_fields_trial = GT.num_fields(trial_space) counter = matrix_strategy.counter(setup) function loop(glue_test,glue_trial,field_per_dim,measure,contribution) + field_test,field_trial = field_per_dim sface_to_faces_test, _,sface_to_faces_around_test = target_face(glue_test) sface_to_faces_trial, _,sface_to_faces_around_trial = target_face(glue_trial) - face_to_dofs_test = face_dofs(test_space,field_per_dim[1]) - face_to_dofs_trial = face_dofs(trial_space,field_per_dim[2]) - field_test,field_trial = field_per_dim + face_to_dofs_test = face_dofs(test_space,field_test) + face_to_dofs_trial = face_dofs(trial_space,field_trial) term_qty = GT.term(contribution) term_npoints = GT.num_points(measure) sface = :sface @@ -415,10 +427,16 @@ function assemble_matrix_fill(state) faces_around_test = sface_to_faces_around_test[sface] faces_around_trial = sface_to_faces_around_trial[sface] for ($face_around_test,face_test) in zip(faces_around_test,faces_test) + if face_test == 0 + continue + end $(s_qty[3]) dofs_test = face_to_dofs_test[face_test] ndofs_test = length(dofs_test) for ($face_around_trial,face_trial) in zip(faces_around_trial,faces_trial) + if face_trial == 0 + continue + end $(s_qty[4]) dofs_trial = face_to_dofs_trial[face_trial] ndofs_trial = length(dofs_trial) diff --git a/src/interpolation.jl b/src/interpolation.jl index 33a94a63..b66bcedf 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -575,6 +575,9 @@ end # TODO a better name function face_shape_function(rid_to_fs,face_to_rid,face,dof) + if dof == 0 || face == 0 + return first(first(rid_to_fs)) + end fs = reference_value(rid_to_fs,face_to_rid,face) fs[dof] end @@ -602,7 +605,7 @@ function shape_functions(a::AbstractSpace,dim) else @term begin fs = face_shape_function($refid_to_funs_sym,$face_to_refid_sym,$face,$dof) - bool = $face_around == $(face_around_per_dim[dim]) + bool = ($face_around == $(face_around_per_dim[dim])) & ($face != 0) mask_function(fs,bool) end end @@ -807,10 +810,17 @@ function shape_functions(a::CartesianProductSpace,dim,field) GT.quantity(GT.prototype(qty),GT.domain(qty)) do index # TODO field_per_dim[dim] has possibly already a numeric value field_per_dim = index.field_per_dim + dof_per_dim = index.dof_per_dim + dof = @term begin + coeff = $field == $(field_per_dim[dim]) + coeff*$(dof_per_dim[dim]) + end + dof_per_dim2 = Base.setindex(dof_per_dim,dof,dim) + index2 = replace_dof_per_dim(index,dof_per_dim2) @assert field_per_dim !== nothing @term begin bool = $field == $(field_per_dim[dim]) - mask_function($(t(index)),bool) + mask_function($(t(index2)),bool) end end end From 7e50045b8ac19f5f2eb8873a5973c258c941b508 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Thu, 29 Aug 2024 16:10:22 +0200 Subject: [PATCH 20/20] Activating examples --- GalerkinToolkitExamples/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GalerkinToolkitExamples/test/runtests.jl b/GalerkinToolkitExamples/test/runtests.jl index aecba88b..8eefcb73 100644 --- a/GalerkinToolkitExamples/test/runtests.jl +++ b/GalerkinToolkitExamples/test/runtests.jl @@ -3,7 +3,7 @@ module GalerkinToolkitExamplesTests using Test @testset "GalerkinToolkitExamples" begin - #@testset "poisson" begin include("poisson_tests.jl") end + @testset "poisson" begin include("poisson_tests.jl") end end end # module