diff --git a/GalerkinToolkitExamples/Project.toml b/GalerkinToolkitExamples/Project.toml index 9a58e0a..8ab5984 100644 --- a/GalerkinToolkitExamples/Project.toml +++ b/GalerkinToolkitExamples/Project.toml @@ -4,6 +4,8 @@ authors = ["Francesc Verdugo "] version = "0.1.0" [deps] +AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GalerkinToolkit = "5e3ba9c4-dd81-444d-b69a-0e7bd7bf60a4" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" @@ -20,14 +22,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" + +[compat] +PartitionedSolvers = "0.3" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test"] - -[compat] -PartitionedSolvers = "0.3" - diff --git a/GalerkinToolkitExamples/src/poisson.jl b/GalerkinToolkitExamples/src/poisson.jl index 0b41542..8d998f3 100644 --- a/GalerkinToolkitExamples/src/poisson.jl +++ b/GalerkinToolkitExamples/src/poisson.jl @@ -13,6 +13,8 @@ using SparseArrays using PartitionedArrays using TimerOutputs using WriteVTK +using AbstractTrees + function main(params_in) params_default = default_params() @@ -26,18 +28,15 @@ function main(params_in) end end -gradient(u) = x->ForwardDiff.gradient(u,x) -jacobian(u) = x->ForwardDiff.jacobian(u,x) laplacian(u,x) = tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x)) -laplacian(u) = x->laplacian(u,x) +# TODO hide this +laplacian(u::GT.AbstractQuantity,x::GT.AbstractQuantity) = GT.call(laplacian,u,x) -Δ(u) = GT.call(laplacian,u) -∇(u) = GT.call(gradient,u) -Δ(u,x) = Δ(u)(x) -∇(u,x) = ∇(u)(x) +const ∇ = ForwardDiff.gradient +const Δ = laplacian -mean(u,x) = 0.5*(u(x)[+]+u(x)[-]) -jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x) +mean(u) = 0.5*(u[1]+u[2]) +jump(u,n) = u[2]*n[2] + u[1]*n[1] function main_automatic(params) timer = params[:timer] @@ -54,7 +53,7 @@ function main_automatic(params) u = GT.analytical_field(params[:u],Ω) f(x) = -Δ(u,x) - n_Γn = GT.unit_normal(Γn,Ω) + n_Γn = GT.unit_normal(mesh,D-1) g(x) = n_Γn(x)⋅∇(u,x) interpolation_degree = params[:interpolation_degree] @@ -68,7 +67,7 @@ function main_automatic(params) GT.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__") Λ = GT.skeleton(mesh;physical_names=["__INTERIOR_FACES__"]) dΛ = GT.measure(Λ,integration_degree) - n_Λ = GT.unit_normal(Λ,Ω) + n_Λ = GT.unit_normal(mesh,D-1) h_Λ = GT.face_diameter_field(Λ) else conformity = :default @@ -79,7 +78,7 @@ function main_automatic(params) uhd = GT.dirichlet_field(Float64,V) GT.interpolate_dirichlet!(u,uhd) else - n_Γd = GT.unit_normal(Γd,Ω) + n_Γd = GT.unit_normal(mesh,D-1) h_Γd = GT.face_diameter_field(Γd) dΓd = GT.measure(Γd,integration_degree) V = GT.lagrange_space(Ω,interpolation_degree;conformity) @@ -98,7 +97,7 @@ function main_automatic(params) end if params[:discretization_method] === :interior_penalty r += ∫( x-> - (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-jump(v,n_Λ,x)⋅mean(∇(u),x)-mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ) + (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-jump(v(x),n_Λ(x))⋅mean(∇(u,x))-mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ) end r end @@ -202,7 +201,7 @@ function main_hand_written(params) dΩ = GT.measure(Ω,integration_degree) u = params[:u] - f(x) = -laplacian(u,x) + f(x) = -Δ(u,x) # Set interpolation space -> Lagrange polynomial basis function, the element type and geometry V = GT.lagrange_space(Ω,interpolation_degree;dirichlet_boundary=Γd) @@ -564,8 +563,8 @@ function default_params() params[:dirichlet_tags] = ["boundary"] params[:neumann_tags] = String[] params[:dirichlet_method] = :strong - params[:integration_degree] = 1 - params[:interpolation_degree] = 2*params[:integration_degree] + params[:interpolation_degree] = 1 + params[:integration_degree] = 2*params[:interpolation_degree] params[:solver] = PS.LinearAlgebra_lu params[:example_path] = joinpath(mkpath(joinpath(@__DIR__,"..","output")),"example_001") params[:export_vtu] = true diff --git a/GalerkinToolkitExamples/test/poisson_tests.jl b/GalerkinToolkitExamples/test/poisson_tests.jl index 92c7b45..f8c02a9 100644 --- a/GalerkinToolkitExamples/test/poisson_tests.jl +++ b/GalerkinToolkitExamples/test/poisson_tests.jl @@ -19,7 +19,7 @@ for discretization_method in (:continuous_galerkin,) params[:dirichlet_tags] = ["1-face-1","1-face-2","1-face-3","1-face-4"] params[:discretization_method] = discretization_method params[:dirichlet_method] = dirichlet_method - params[:integration_degree] = 2 + params[:interpolation_degree] = 2 results = Poisson.main(params) @test results[:error_l2_norm] < tol @test results[:error_h1_norm] < tol @@ -35,7 +35,7 @@ for discretization_method in (:interior_penalty,:continuous_galerkin) params[:neumann_tags] = ["1-face-2"] params[:discretization_method] = discretization_method params[:dirichlet_method] = dirichlet_method - params[:integration_degree] = 2 + params[:interpolation_degree] = 2 results = Poisson.main(params) @test results[:error_l2_norm] < tol @test results[:error_h1_norm] < tol diff --git a/docs/src/examples/methods.jl b/docs/src/examples/methods.jl index a53a8da..b9aa69f 100644 --- a/docs/src/examples/methods.jl +++ b/docs/src/examples/methods.jl @@ -11,21 +11,22 @@ import FileIO # hide domain = (0,1,0,1) cells = (10,10) mesh = GT.cartesian_mesh(domain,cells) +D = GT.num_dims(mesh) GT.label_boundary_faces!(mesh;physical_name="Γd") Ω = GT.interior(mesh) Γd = GT.boundary(mesh;physical_names=["Γd"]) -n_Γd = GT.unit_normal(Γd,Ω) +n_Γd = GT.unit_normal(mesh,D-1) u = GT.analytical_field(x->sum(x),Ω) #TODO these two definitions should be enough -#∇ = ForwardDiff.gradient -#Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x)) -gradient(u) = x->ForwardDiff.gradient(u,x) -jacobian(u) = x->ForwardDiff.jacobian(u,x) -laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x)) -Δ(u) = GT.call(laplacian,u) -∇(u) = GT.call(gradient,u) -Δ(u,x) = Δ(u)(x) -∇(u,x) = ∇(u)(x) +∇ = ForwardDiff.gradient +Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x)) +#gradient(u) = x->ForwardDiff.gradient(u,x) +#jacobian(u) = x->ForwardDiff.jacobian(u,x) +#laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x)) +#Δ(u) = GT.call(laplacian,u) +#∇(u) = GT.call(gradient,u) +#Δ(u,x) = Δ(u)(x) +#∇(u,x) = ∇(u)(x) tol = 1e-10 order = 2 γ = order*(order+1)/10 @@ -54,7 +55,7 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ))) # ## Nitche Method V = GT.lagrange_space(Ω,order) -n_Γd = GT.unit_normal(Γd,Ω) +n_Γd = GT.unit_normal(mesh,D-1) h_Γd = GT.face_diameter_field(Γd) dΓd = GT.measure(Γd,2*order) a_Γd = (u,v) -> GT.∫( @@ -90,18 +91,18 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ))) # ## Interior Penalty -mean(u,x) = 0.5*(u(x)[+]+u(x)[-]) -jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x) +mean(u) = 0.5*(u[1]+u[2]) +jump(u,n) = u[2]*n[2] + u[1]*n[1] GT.label_interior_faces!(mesh;physical_name="Λ") Λ = GT.skeleton(mesh;physical_names=["Λ"]) dΛ = GT.measure(Λ,2*order) -n_Λ = GT.unit_normal(Λ,Ω) +n_Λ = GT.unit_normal(mesh,D-1) h_Λ = GT.face_diameter_field(Λ) V = GT.lagrange_space(Ω,order;conformity=:L2) a_Λ = (u,v) -> GT.∫( - x-> (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)- - jump(v,n_Λ,x)⋅mean(∇(u),x)- - mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ) + x-> (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))- + jump(v(x),n_Λ(x))⋅mean(∇(u,x))- + mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ) a_ip = (u,v) -> a_Ω(u,v) + a_Γd(u,v) + a_Λ(u,v) l_ip = l_ni p = GT.linear_problem(Float64,V,a_ip,l_ip) diff --git a/src/assembly.jl b/src/assembly.jl index 5dd1d63..7cd8f51 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -192,11 +192,17 @@ function assemble_vector_count(integral,state) for field in 1:nfields Dface_to_sDface = field_to_Dface_to_sDface[field] D = field_to_D[field] + if D < d + continue + end face_to_Dfaces = face_incidence(topo,d,D) sDface_to_dofs = field_to_sDface_to_dofs[field] Dfaces = face_to_Dfaces[face] for Dface in Dfaces sDface = Dface_to_sDface[Dface] + if sDface == 0 + continue + end dofs = sDface_to_dofs[sDface] for dof in dofs counter = vector_strategy.count(counter,setup,dof,field) @@ -230,7 +236,7 @@ function assemble_vector_fill!(integral,state) nfields = GT.num_fields(space) field_to_domain = map(field->domain(space,field),1:nfields) field_to_sDface_to_dofs = map(field->face_dofs(space,field),1:nfields) - field_to_Dface_to_sDface = map(faces,field_to_domain) + field_to_Dface_to_sDface = map(inverse_faces,field_to_domain) field_to_D = map(num_dims,field_to_domain) counter = vector_strategy.counter(setup) T = vector_strategy.scalar_type(setup) @@ -270,6 +276,9 @@ function assemble_vector_fill!(integral,state) $(s_qty[3]) Dface_to_sDface = args.field_to_Dface_to_sDface[$field] D = args.field_to_D[$field] + if D < args.d + continue + end face_to_Dfaces = face_incidence(args.topo,args.d,D) sDface_to_dofs = args.field_to_sDface_to_dofs[$field] Dfaces = face_to_Dfaces[$face] @@ -277,6 +286,9 @@ function assemble_vector_fill!(integral,state) $(s_qty[4]) fill!(b,zero(eltype(b))) sDface = Dface_to_sDface[Dface] + if sDface == 0 + continue + end dofs = sDface_to_dofs[sDface] for $point in 1:npoints $(s_qty[5]) @@ -371,7 +383,7 @@ function assemble_matrix_count(integral,state) axis_to_nfields = map(num_fields,axis_to_space) axis_to_field_to_domain = map(space->map(field->domain(space,field),1:num_fields(space)),axis_to_space) axis_to_field_to_sDface_to_dofs = map(space->map(field->face_dofs(space,field),1:num_fields(space)),axis_to_space) - axis_to_field_to_Dface_to_sDface = map(field_to_domain->map(faces,field_to_domain),axis_to_field_to_domain) + axis_to_field_to_Dface_to_sDface = map(field_to_domain->map(inverse_faces,field_to_domain),axis_to_field_to_domain) # TODO: faces or inverse_faces?? axis_to_field_to_D = map(field_to_domain->map(num_dims,field_to_domain),axis_to_field_to_domain) counter = matrix_strategy.counter(setup) for measure_and_contribution in contributions @@ -384,26 +396,38 @@ function assemble_matrix_count(integral,state) for sface in 1:nsfaces face = sface_to_face[sface] for field_test in 1:axis_to_nfields[test] - test_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[field_test][test] - test_D = axis_to_field_to_D[field_test][test] + test_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[test][field_test] + test_D = axis_to_field_to_D[test][field_test] + if test_D < d + continue + end test_face_to_Dfaces = face_incidence(topo,d,test_D) - test_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[field_test][test] + test_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[test][field_test] test_Dfaces = test_face_to_Dfaces[face] for field_trial in 1:axis_to_nfields[trial] - trial_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[field_trial][trial] - trial_D = axis_to_field_to_D[field_trial][trial] + trial_Dface_to_sDface = axis_to_field_to_Dface_to_sDface[trial][field_trial] + trial_D = axis_to_field_to_D[trial][field_trial] + if trial_D < d + continue + end trial_face_to_Dfaces = face_incidence(topo,d,trial_D) - trial_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[field_trial][trial] + trial_sDface_to_dofs = axis_to_field_to_sDface_to_dofs[trial][field_trial] trial_Dfaces = trial_face_to_Dfaces[face] - for trial_Dface in trial_Dfaces #TODO Typo this should be test - trial_sDface = trial_Dface_to_sDface[trial_Dface] - trial_dofs = trial_sDface_to_dofs[trial_sDface] + for test_Dface in test_Dfaces + test_sDface = test_Dface_to_sDface[test_Dface] + if test_sDface == 0 + continue + end + test_dofs = test_sDface_to_dofs[test_sDface] for trial_Dface in trial_Dfaces trial_sDface = trial_Dface_to_sDface[trial_Dface] + if trial_sDface == 0 + continue + end trial_dofs = trial_sDface_to_dofs[trial_sDface] for test_dof in test_dofs for trial_dof in trial_dofs - counter = matrix_strategy.count(counter,setup,test_dof,trial_dof,test_field,trial_field) + counter = matrix_strategy.count(counter,setup,test_dof,trial_dof,field_test,field_trial) end end end @@ -422,106 +446,116 @@ function assemble_matrix_allocate(state) end function assemble_matrix_fill!(integral,state) - (;test_space,trial_space,fields_test,fields_trial,alloc,glues_test,glues_trial,matrix_strategy,setup) = state + (;test_space,trial_space,alloc,matrix_strategy,setup) = state + contributions = GT.contributions(integral) - num_fields_test = GT.num_fields(test_space) - num_fields_trial = GT.num_fields(trial_space) + test, trial = 1, 2 + axis_to_space = (test_space,trial_space) + axis_to_nfields = map(num_fields,axis_to_space) + axis_to_field_to_domain = map(space->map(field->domain(space,field),1:num_fields(space)),axis_to_space) + axis_to_field_to_sDface_to_dofs = map(space->map(field->face_dofs(space,field),1:num_fields(space)),axis_to_space) + axis_to_field_to_Dface_to_sDface = map(field_to_domain->map(inverse_faces,field_to_domain),axis_to_field_to_domain) # TODO: faces or inverse_faces?? + axis_to_field_to_D = map(field_to_domain->map(num_dims,field_to_domain),axis_to_field_to_domain) 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_test) - face_to_dofs_trial = face_dofs(trial_space,field_trial) - term_qty = GT.term(contribution) + T = matrix_strategy.scalar_type(setup) + b = zeros(T,max_local_dofs(trial_space), max_local_dofs(test_space)) + + for measure_and_contribution in contributions + measure,contribution = measure_and_contribution + domain = GT.domain(measure) + sface_to_face = faces(domain) + topo = topology(mesh(domain)) + d = num_dims(domain) + form_arity = 2 + index = GT.generate_index(domain,form_arity) + t = GT.term(contribution,index) 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_qty = t |> expression |> simplify + expr_npoints = term_npoints(index) |> expression |> simplify + face = face_index(index,d) + point = point_index(index) + axes = (test, trial) + + idof_test, idof_trial = map(axis -> dof_index(index,axis), axes) # TODO: test, trial + test_Dface_around, trial_Dface_around = map(axis -> face_around_index(index,axis), axes) + field_test, field_trial = map(axis -> field_index(index,axis), axes) + + s_qty = GT.topological_sort(expr_qty,(face, field_test, field_trial, test_Dface_around, trial_Dface_around, point, idof_test, idof_trial)) # 8 deps + s_npoints = GT.topological_sort(expr_npoints,(face,)) + expr = quote 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,field_test,field_trial, - alloc,matrix_strategy,setup) = args - $(unpack_storage(index.dict,:storage)) + $(unpack_index_storage(index,:storage)) $(s_qty[1]) $(s_npoints[1]) - nsfaces = length(sface_to_faces_test) - z = zero(matrix_strategy.scalar_type(setup)) - for $sface in 1:nsfaces + b = args.b + nsfaces = length(args.sface_to_face) + for sface in 1:nsfaces + $face = args.sface_to_face[sface] $(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) - if face_test == 0 + for $field_test in 1:args.axis_to_nfields[$test] + $(s_qty[3]) + test_Dface_to_sDface = args.axis_to_field_to_Dface_to_sDface[$test][$field_test] + test_D = args.axis_to_field_to_D[$test][$field_test] + if test_D < args.d 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 + test_face_to_Dfaces = face_incidence(args.topo,args.d,test_D) + test_sDface_to_dofs = args.axis_to_field_to_sDface_to_dofs[$test][$field_test] + test_Dfaces = test_face_to_Dfaces[$face] + for $field_trial in 1:args.axis_to_nfields[$trial] + $(s_qty[4]) + trial_Dface_to_sDface = args.axis_to_field_to_Dface_to_sDface[$trial][$field_trial] + trial_D = args.axis_to_field_to_D[$trial][$field_trial] + if trial_D < args.d continue end - $(s_qty[4]) - dofs_trial = face_to_dofs_trial[face_trial] - ndofs_trial = length(dofs_trial) - for $idof_test in 1:ndofs_test + trial_face_to_Dfaces = face_incidence(args.topo,args.d,trial_D) + trial_sDface_to_dofs = args.axis_to_field_to_sDface_to_dofs[$trial][$field_trial] + trial_Dfaces = trial_face_to_Dfaces[$face] + for ($test_Dface_around,test_Dface) in enumerate(test_Dfaces) $(s_qty[5]) - dof_test = dofs_test[$idof_test] - for $idof_trial in 1:ndofs_trial + test_sDface = test_Dface_to_sDface[test_Dface] + if test_sDface == 0 + continue + end + test_dofs = test_sDface_to_dofs[test_sDface] + for ($trial_Dface_around,trial_Dface) in enumerate(trial_Dfaces) $(s_qty[6]) - dof_trial = dofs_trial[$idof_trial] - v = z + trial_sDface = trial_Dface_to_sDface[trial_Dface] + if trial_sDface == 0 + continue + end + trial_dofs = trial_sDface_to_dofs[trial_sDface] + fill!(b,zero(eltype(b))) for $point in 1:npoints - v += $(s_qty[7]) + $(s_qty[7]) + for ($idof_test, dof_test) in enumerate(test_dofs) + $(s_qty[8]) + for ($idof_trial, dof_trial) in enumerate(trial_dofs) + b[$idof_test, $idof_trial] += $(s_qty[9]) + end + end + end + for ($idof_test, dof_test) in enumerate(test_dofs) + for ($idof_trial, dof_trial) in enumerate(trial_dofs) + counter = args.matrix_strategy.set!(args.alloc,counter,args.setup,b[$idof_test, $idof_trial],dof_test,dof_trial,$field_test,$field_trial) + end end - counter = matrix_strategy.set!(alloc,counter,setup,v,dof_test,dof_trial,field_test,field_trial) end end end end + end 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,field_test,field_trial, - alloc,matrix_strategy,setup) + storage = index_storage(index) + args = (;d,axis_to_nfields,alloc,matrix_strategy,setup,b,topo,sface_to_face,axis_to_field_to_D,axis_to_field_to_Dface_to_sDface,axis_to_field_to_sDface_to_dofs) 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) - 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 - glue_trial = field_to_glue_trial[field_trial] - if glue_test === nothing || glue_trial == nothing - continue - end - field_per_dim = (field_test,field_trial) - loop(glue_test,glue_trial,field_per_dim,measure,contribution) - end - end end state end diff --git a/src/interpolation.jl b/src/interpolation.jl index 6874cc5..cf912d4 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -1004,7 +1004,7 @@ function form_argument(a::AbstractSpace,axis,field=1) refid_to_reffes = reference_fes(a) refid_to_funs = map(GT.shape_functions,refid_to_reffes) domain = GT.domain(a) - qty = form_argument(axis,field,refid_to_funs,domain;reference=true,face_reference_id=face_to_refid) + qty = form_argument(axis,field,refid_to_funs,domain;reference=true)# TODO,face_reference_id=face_to_refid) if is_reference_domain(GT.domain(a)) return qty end diff --git a/src/symbolics.jl b/src/symbolics.jl index 2a97019..34967b2 100644 --- a/src/symbolics.jl +++ b/src/symbolics.jl @@ -332,6 +332,15 @@ function binary_call_term(f,a,b::BoundaryTerm) boundary_term(b.dim,b.dim_around,t,b.face_around) end +function binary_call_term(f,a::BoundaryTerm,b) + t = binary_call_term(f,a.term,b) + boundary_term(a.dim,a.dim_around,t,a.face_around) +end + +function binary_call_term(f,a::BoundaryTerm,b::BoundaryTerm) + BinaryCallTerm(f,a,b) +end + skeleton_term(args...) = SkeletonTerm(args...) function skeleton_term(dim,dim_around,term::BoundaryTerm) @@ -825,7 +834,10 @@ end free_dims(a::PhysicalMapTerm) = [a.dim] -prototype(a::PhysicalMapTerm) = x-> zero(SVector{val_parameter(a.dim),Float64}) +function prototype(a::PhysicalMapTerm) + D = num_ambient_dims(mesh(index(a))) + x-> zero(SVector{D,Float64}) +end function Base.:(==)(a::PhysicalMapTerm,b::PhysicalMapTerm) @@ -964,6 +976,28 @@ function binary_call_term(d::typeof(ForwardDiff.gradient),a::ComposedWithInverse call(\,Jt,gradf) end +function binary_call_term(d::typeof(ForwardDiff.jacobian),a::ComposedWithInverseTerm,b::FunctionCallTerm) + phi1 = a.arg2.arg1 + phi2 = b.arg1 + if phi1 != phi2 + return binary_call_term_physical_maps(d,a,b,phi1,phi2) + end + phi = phi1 + f = a.arg1 + x = b.arg2 + J = call(ForwardDiff.jacobian,phi,x) + Jf = call(d,f,x) + call(/,Jf,J) +end + +# w(q) = u(ϕ(q)) +# ∇(w,q) = Jt(ϕ,q)*∇(u,ϕ(q)) +# Jt(ϕ,q)\∇(w,q) = ∇(u,ϕ(q)) +# +# +# J(w,q) = J(u,ϕ(q))*J(ϕ,q) +# J(w,q)/J(ϕ,q) = J(u,ϕ(q)) + function binary_call_term_physical_maps(op,a,b,phi1,phi2) return BinaryCallTerm(op,a,b) end @@ -976,7 +1010,7 @@ function binary_call_term_physical_maps(::typeof(call),a,b,phi1::PhysicalMapTerm f(phi(x)) end -function binary_call_term_physical_maps(::typeof(ForwardDiff.gradient),a,b,phi1::PhysicalMapTerm,phi2::PhysicalMapTerm) +function binary_call_term_physical_maps(d::typeof(ForwardDiff.gradient),a,b,phi1::PhysicalMapTerm,phi2::PhysicalMapTerm) @assert phi1.dim != phi2.dim phi = reference_map_term(phi2.dim,phi1.dim,index(a)) f = a.arg1 @@ -1237,7 +1271,7 @@ struct UnitNormalTerm{A} <: AbstractTerm end free_dims(a::UnitNormalTerm) = free_dims(a.outwards_normals) -prototype(a::UnitNormalTerm) = prototype(a.outwards_normals) +prototype(a::UnitNormalTerm) = x -> prototype(a.outwards_normals) index(a::UnitNormalTerm) = index(a.outwards_normals) function expression(c::BinaryCallTerm{typeof(call),<:UnitNormalTerm,<:Any}) diff --git a/test/assembly_tests.jl b/test/assembly_tests.jl index a345199..42a5c8a 100644 --- a/test/assembly_tests.jl +++ b/test/assembly_tests.jl @@ -289,9 +289,8 @@ mesh = GT.cartesian_mesh(domain,cells) k = 1 V = GT.lagrange_space(Ω,k) dΩ = GT.measure(Ω,2*k) -gradient(u) = x->ForwardDiff.gradient(u,x) -∇(u,x) = GT.call(gradient,u)(x) -a(u,v) = GT.∫( x->∇(u,x)⋅∇(v,x), dΩ) +∇3 = ForwardDiff.gradient +a(u,v) = GT.∫( x->∇3(u,x)⋅∇3(v,x), dΩ) l(v) = 0 p = GT.linear_problem(Float64,V,a,l) PS.matrix(p) |> display