From 62bac62ba6dd7e34493781556cefdf4c0800b32f Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Mon, 5 Feb 2024 18:30:45 +0100 Subject: [PATCH] Cleanup in example001 --- .../GalerkinToolkitExamples/src/example001.jl | 444 ++++++++++++++---- .../test/example001_tests.jl | 9 +- 2 files changed, 373 insertions(+), 80 deletions(-) diff --git a/extensions/GalerkinToolkitExamples/src/example001.jl b/extensions/GalerkinToolkitExamples/src/example001.jl index 957a5256..03119a85 100644 --- a/extensions/GalerkinToolkitExamples/src/example001.jl +++ b/extensions/GalerkinToolkitExamples/src/example001.jl @@ -6,6 +6,7 @@ using StaticArrays using LinearAlgebra using SparseArrays using WriteVTK +using PartitionedArrays: JaggedArray, dense_vector # This one implements a vanilla sequential iso-parametric Poisson solver by only # using the mesh interface. @@ -118,7 +119,7 @@ function setup_isomap(params,integration) end rid_to_shape_vals = map(first,shape_funs) rid_to_shape_grads = map(last,shape_funs) - face_to_nodes = gk.face_nodes(mesh,d) + face_to_nodes = JaggedArray(gk.face_nodes(mesh,d)) node_to_coords = gk.node_coordinates(mesh) isomap = (;face_to_nodes,node_to_coords,face_to_rid,rid_to_shape_vals,rid_to_shape_grads,d) end @@ -142,6 +143,28 @@ function setup_user_funs(params) user_funs = (;u,f,g) end +function setup_dofs(params,dirichlet_bcs,cell_isomap,face_isomap) + free_and_dirichlet_nodes = dirichlet_bcs.free_and_dirichlet_nodes + node_to_free_node = gk.permutation(free_and_dirichlet_nodes) + n_free = length(first(free_and_dirichlet_nodes)) + n_dofs = n_free + cell_to_nodes = cell_isomap.face_to_nodes + face_to_nodes = face_isomap.face_to_nodes + function node_to_dof(node) + free_node = node_to_free_node[node] + if free_node <= n_free + return Int(free_node) + end + Int(-(free_node-n_free)) + end + cell_to_dofs_data = node_to_dof.(cell_to_nodes.data) + face_to_dofs_data = node_to_dof.(face_to_nodes.data) + cell_to_dofs = JaggedArray(cell_to_dofs_data,cell_to_nodes.ptrs) + face_to_dofs = JaggedArray(face_to_dofs_data,face_to_nodes.ptrs) + dofs = (;cell_to_dofs,face_to_dofs,n_dofs) + dofs +end + function setup(params) dirichlet_bcs = setup_dirichlet_bcs(params) neumann_bcs = setup_neumann_bcs(params) @@ -149,9 +172,10 @@ function setup(params) face_integration = setup_integration(params,:faces) cell_isomap = setup_isomap(params,cell_integration) face_isomap = setup_isomap(params,face_integration) + dofs = setup_dofs(params,dirichlet_bcs,cell_isomap,face_isomap) user_funs = setup_user_funs(params) solver = params[:solver] - state = (;solver,dirichlet_bcs,neumann_bcs,cell_integration,cell_isomap,face_integration,face_isomap,user_funs) + state = (;solver,dirichlet_bcs,neumann_bcs,cell_integration,cell_isomap,face_integration,face_isomap,user_funs,dofs) end function add_basic_info(results,params,state) @@ -167,7 +191,79 @@ function add_basic_info(results,params,state) results end -function assemble_system_loop(state) +function assemble_system(state) + args = assemble_sytem_coo(state) + I,J,V,II,VV = args + n_dofs = state.dofs.n_dofs + A = sparse(I,J,V,n_dofs,n_dofs) + b = dense_vector(II,VV,n_dofs) + A,b +end + +function assemble_sytem_coo(state) + args = assemble_sybmolic(state) + ii_coo = assemble_numeric_cells!(args...,state) + assemble_numeric_faces!(args...,state,ii_coo) + set_dirichlet_bcs!(args...,state) + filter_dirichlet!(args...) + args +end + +function assemble_sybmolic(state) + cell_to_dofs = state.dofs.cell_to_dofs + face_to_dofs = state.dofs.face_to_dofs + neum_face_to_face = state.neumann_bcs.neum_face_to_face + + n_coo = 0 + nn_coo = 0 + ncells = length(cell_to_dofs) + nfaces = length(face_to_dofs) + for cell in 1:ncells + dofs = cell_to_dofs[cell] + ndofs = length(dofs) + n_coo += ndofs*ndofs + nn_coo += ndofs + end + for face in neum_face_to_face + dofs = face_to_dofs[face] + ndofs = length(dofs) + nn_coo += ndofs + end + + I_coo = zeros(Int,n_coo) + J_coo = zeros(Int,n_coo) + V_coo = zeros(Float64,n_coo) + II_coo = zeros(Int,nn_coo) + VV_coo = zeros(Float64,nn_coo) + + n_coo = 0 + nn_coo = 0 + for cell in 1:ncells + dofs = cell_to_dofs[cell] + ndofs = length(dofs) + for i in 1:ndofs + nn_coo += 1 + II_coo[nn_coo] = dofs[i] + for j in 1:ndofs + n_coo += 1 + I_coo[n_coo] = dofs[i] + J_coo[n_coo] = dofs[j] + end + end + end + for face in neum_face_to_face + dofs = face_to_dofs[face] + ndofs = length(dofs) + for i in 1:ndofs + nn_coo += 1 + II_coo[nn_coo] = dofs[i] + end + end + + I_coo,J_coo,V_coo,II_coo,VV_coo +end + +function assemble_numeric_cells!(I_coo,J_coo,V_coo,II_coo,VV_coo,state) cell_to_nodes = state.cell_isomap.face_to_nodes cell_to_rid = state.cell_isomap.face_to_rid @@ -179,49 +275,22 @@ function assemble_system_loop(state) f = state.user_funs.f w = state.cell_integration.rid_to_weights d = state.cell_integration.d - - # Count coo entries - n_coo = 0 ncells = length(cell_to_nodes) - node_to_free_node = gk.permutation(free_and_dirichlet_nodes) - nfree = length(first(free_and_dirichlet_nodes)) - for cell in 1:ncells - nodes = cell_to_nodes[cell] - for nj in nodes - if node_to_free_node[nj] > nfree - continue - end - for ni in nodes - if node_to_free_node[ni] > nfree - continue - end - n_coo += 1 - end - end - end - - # Allocate coo values - I_coo = zeros(Int32,n_coo) - J_coo = zeros(Int32,n_coo) - V_coo = zeros(Float64,n_coo) - b = zeros(nfree) # Allocate auxiliary buffers ∇x = map(i->similar(i,size(i,2)),∇s) Aes = map(i->zeros(size(i,2),size(i,2)),∇s) fes = map(i->zeros(size(i,2)),∇s) - ues = map(i->zeros(size(i,2)),∇s) Tx = SVector{d,Float64} TJ = typeof(zero(Tx)*zero(Tx)') - # Fill coo values i_coo = 0 + ii_coo = 0 for cell in 1:ncells rid = cell_to_rid[cell] nodes = cell_to_nodes[cell] Ae = Aes[rid] - ue = ues[rid] fe = fes[rid] nl = length(nodes) ∇se = ∇s[rid] @@ -231,15 +300,7 @@ function assemble_system_loop(state) nq = length(we) fill!(Ae,zero(eltype(Ae))) fill!(fe,zero(eltype(Ae))) - fill!(ue,zero(eltype(Ae))) - for k in 1:nl - nk = nodes[k] - gk = node_to_free_node[nk] - if gk <= nfree - continue - end - ue[k] = u_dirichlet[gk-nfree] - end + # Integrate for iq in 1:nq Jt = zero(TJ) xint = zero(Tx) @@ -266,54 +327,34 @@ function assemble_system_loop(state) fe[k] += fx*se[iq,k]*dV end end + # Set the result in the output array for i in 1:nl - for j in 1:nl - fe[i] -= Ae[i,j]*ue[j] - end - end - for i in 1:nl - ni = nodes[i] - gi = node_to_free_node[ni] - if gi > nfree - continue - end - b[gi] += fe[i] + ii_coo += 1 + VV_coo[ii_coo] = fe[i] end for j in 1:nl - nj = nodes[j] - gj = node_to_free_node[nj] - if gj > nfree - continue - end for i in 1:nl - ni = nodes[i] - gi = node_to_free_node[ni] - if gi > nfree - continue - end i_coo += 1 - I_coo[i_coo] = gi - J_coo[i_coo] = gj V_coo[i_coo] = Ae[i,j] end end end - add_neumann_bcs!(b,state) - I_coo,J_coo,V_coo,b + + ii_coo end -function add_neumann_bcs!(b,state) +function assemble_numeric_faces!(I_coo,J_coo,V_coo,II_coo,VV_coo,state,ii_coo) neum_face_to_face = state.neumann_bcs.neum_face_to_face if length(neum_face_to_face) == 0 return nothing end + face_to_rid = state.face_isomap.face_to_rid face_to_nodes = state.face_isomap.face_to_nodes ∇s_f = state.face_isomap.rid_to_shape_grads s_f = state.face_isomap.rid_to_shape_vals node_to_x = state.face_isomap.node_to_coords - free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes g = state.user_funs.g w_f = state.face_integration.rid_to_weights d = state.cell_integration.d @@ -321,8 +362,6 @@ function add_neumann_bcs!(b,state) fes = map(i->zeros(size(i,2)),s_f) Tx = SVector{d,Float64} Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)') - node_to_free_node = gk.permutation(free_and_dirichlet_nodes) - nfree = length(first(free_and_dirichlet_nodes)) for face in neum_face_to_face rid = face_to_rid[face] @@ -352,23 +391,63 @@ function add_neumann_bcs!(b,state) end end for i in 1:nl - ni = nodes[i] - gi = node_to_free_node[ni] - if gi > nfree + ii_coo += 1 + VV_coo[ii_coo] = fe[i] + end + end +end + +function set_dirichlet_bcs!(I_coo,J_coo,V_coo,II_coo,VV_coo,state) + + cell_to_dofs = state.dofs.cell_to_dofs + ncells = length(cell_to_dofs) + u_dirichlet = state.dirichlet_bcs.u_dirichlet + i_coo = 0 + ii_coo = 0 + for cell in 1:ncells + dofs = cell_to_dofs[cell] + ndofs = length(dofs) + for i in 1:ndofs + if !(dofs[i]>0) continue end - b[gi] += fe[i] + for j in 1:ndofs + if dofs[j]>0 + continue + end + uj = u_dirichlet[-dofs[j]] + ij = (j-1)*ndofs + i + Aij = V_coo[ij+i_coo] + VV_coo[i+ii_coo] -= Aij*uj + end end + i_coo += ndofs*ndofs + ii_coo += ndofs end - nothing end -function assemble_system(state) - I,J,V,b = assemble_system_loop(state) - nfree = length(b) - A = sparse(I,J,V,nfree,nfree) - A,b +function filter_dirichlet!(I_coo,J_coo,V_coo,II_coo,VV_coo) + # TODO this is only needed when the functions in charge of + # compressing the matrix/vector are not able to skip negative + # values + # NB. It will not work for empty problems + for p in 1:length(I_coo) + i = I_coo[p] + j = J_coo[p] + if i<=0 || j<=0 + I_coo[p] = 1 + J_coo[p] = 1 + V_coo[p] = 0 + end + end + for p in 1:length(II_coo) + i = II_coo[p] + if i<=0 + II_coo[p] = 1 + VV_coo[p] = 0 + end + end end function solve_system(A,b,params) @@ -481,4 +560,211 @@ function export_results(uh,params,state) end end +### From here, it is old stuff to to be deleted + +function assemble_system_loop(state) + + cell_to_nodes = state.cell_isomap.face_to_nodes + cell_to_rid = state.cell_isomap.face_to_rid + free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes + node_to_x = state.cell_isomap.node_to_coords + ∇s = state.cell_isomap.rid_to_shape_grads + s = state.cell_isomap.rid_to_shape_vals + u_dirichlet = state.dirichlet_bcs.u_dirichlet + f = state.user_funs.f + w = state.cell_integration.rid_to_weights + d = state.cell_integration.d + + # Count coo entries + n_coo = 0 + ncells = length(cell_to_nodes) + node_to_free_node = gk.permutation(free_and_dirichlet_nodes) + nfree = length(first(free_and_dirichlet_nodes)) + for cell in 1:ncells + nodes = cell_to_nodes[cell] + for nj in nodes + if node_to_free_node[nj] > nfree + continue + end + for ni in nodes + if node_to_free_node[ni] > nfree + continue + end + n_coo += 1 + end + end + end + + # Allocate coo values + I_coo = zeros(Int32,n_coo) + J_coo = zeros(Int32,n_coo) + V_coo = zeros(Float64,n_coo) + b = zeros(nfree) + + # Allocate auxiliary buffers + ∇x = map(i->similar(i,size(i,2)),∇s) + Aes = map(i->zeros(size(i,2),size(i,2)),∇s) + fes = map(i->zeros(size(i,2)),∇s) + ues = map(i->zeros(size(i,2)),∇s) + + Tx = SVector{d,Float64} + TJ = typeof(zero(Tx)*zero(Tx)') + + # Fill coo values + i_coo = 0 + for cell in 1:ncells + rid = cell_to_rid[cell] + nodes = cell_to_nodes[cell] + Ae = Aes[rid] + ue = ues[rid] + fe = fes[rid] + nl = length(nodes) + ∇se = ∇s[rid] + ∇xe = ∇x[rid] + se = s[rid] + we = w[rid] + nq = length(we) + fill!(Ae,zero(eltype(Ae))) + fill!(fe,zero(eltype(Ae))) + fill!(ue,zero(eltype(Ae))) + for k in 1:nl + nk = nodes[k] + gk = node_to_free_node[nk] + if gk <= nfree + continue + end + ue[k] = u_dirichlet[gk-nfree] + end + for iq in 1:nq + Jt = zero(TJ) + xint = zero(Tx) + for k in 1:nl + x = node_to_x[nodes[k]] + ∇sqx = ∇se[iq,k] + sqx = se[iq,k] + Jt += ∇sqx*x' + xint += sqx*x + end + detJt = det(Jt) + invJt = inv(Jt) + dV = abs(detJt)*we[iq] + for k in 1:nl + ∇xe[k] = invJt*∇se[iq,k] + end + for j in 1:nl + for i in 1:nl + Ae[i,j] += ∇xe[i]⋅∇xe[j]*dV + end + end + fx = f(xint) + for k in 1:nl + fe[k] += fx*se[iq,k]*dV + end + end + for i in 1:nl + for j in 1:nl + fe[i] -= Ae[i,j]*ue[j] + end + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + b[gi] += fe[i] + end + for j in 1:nl + nj = nodes[j] + gj = node_to_free_node[nj] + if gj > nfree + continue + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + i_coo += 1 + I_coo[i_coo] = gi + J_coo[i_coo] = gj + V_coo[i_coo] = Ae[i,j] + end + end + end + add_neumann_bcs!(b,state) + I_coo,J_coo,V_coo,b +end + +function add_neumann_bcs!(b,state) + + neum_face_to_face = state.neumann_bcs.neum_face_to_face + if length(neum_face_to_face) == 0 + return nothing + end + face_to_rid = state.face_isomap.face_to_rid + face_to_nodes = state.face_isomap.face_to_nodes + ∇s_f = state.face_isomap.rid_to_shape_grads + s_f = state.face_isomap.rid_to_shape_vals + node_to_x = state.face_isomap.node_to_coords + free_and_dirichlet_nodes = state.dirichlet_bcs.free_and_dirichlet_nodes + g = state.user_funs.g + w_f = state.face_integration.rid_to_weights + d = state.cell_integration.d + + fes = map(i->zeros(size(i,2)),s_f) + Tx = SVector{d,Float64} + Ts = typeof(zero(SVector{d-1,Float64})*zero(Tx)') + node_to_free_node = gk.permutation(free_and_dirichlet_nodes) + nfree = length(first(free_and_dirichlet_nodes)) + + for face in neum_face_to_face + rid = face_to_rid[face] + nodes = face_to_nodes[face] + ∇se = ∇s_f[rid] + se = s_f[rid] + we = w_f[rid] + nq = length(we) + fe = fes[rid] + fill!(fe,zero(eltype(fe))) + nl = length(nodes) + for iq in 1:nq + Jt = zero(Ts) + xint = zero(Tx) + for k in 1:nl + x = node_to_x[nodes[k]] + ∇sqx = ∇se[iq,k] + sqx = se[iq,k] + Jt += ∇sqx*x' + xint += sqx*x + end + J = transpose(Jt) + dS = sqrt(det(Jt*J))*we[iq] + gx = g(xint) + for k in 1:nl + fe[k] += gx*se[iq,k]*dS + end + end + for i in 1:nl + ni = nodes[i] + gi = node_to_free_node[ni] + if gi > nfree + continue + end + b[gi] += fe[i] + end + end + + nothing +end + +function assemble_system_old(state) + I,J,V,b = assemble_system_loop(state) + nfree = length(b) + A = sparse(I,J,V,nfree,nfree) + A,b +end + + end # module diff --git a/extensions/GalerkinToolkitExamples/test/example001_tests.jl b/extensions/GalerkinToolkitExamples/test/example001_tests.jl index ceda93e7..5751fb7b 100644 --- a/extensions/GalerkinToolkitExamples/test/example001_tests.jl +++ b/extensions/GalerkinToolkitExamples/test/example001_tests.jl @@ -5,8 +5,15 @@ using GalerkinToolkitExamples: Example001 using Test tol = 1.0e-10 + +params = Dict{Symbol,Any}() +params[:mesh] = gk.cartesian_mesh((0,1,0,1),(3,3)) +results = Example001.main(params) +@test results[:eh1] < tol +@test results[:el2] < tol + params = Dict{Symbol,Any}() -params[:mesh] = gk.cartesian_mesh((0,3,0,2),(10,5),complexify=false) +params[:mesh] = gk.cartesian_mesh((0,3,0,2),(5,10),complexify=false) params[:dirichlet_tags] = ["1-face-1","1-face-3","1-face-4"] params[:neumann_tags] = ["1-face-2"] params[:u] = (x) -> sum(x)