diff --git a/.gitignore b/.gitignore index 53aeca3..e8ab185 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,6 @@ Manifest.toml *.e.*.* *.out build/ +*.exo +*.exo.* +*.msh diff --git a/examples/cooks_membrane/generate_meshes.sh b/examples/cooks_membrane/generate_meshes.sh new file mode 100755 index 0000000..8b145b8 --- /dev/null +++ b/examples/cooks_membrane/generate_meshes.sh @@ -0,0 +1,2 @@ +gmsh geometry.geo -2 -setnumber element_order 1 -o mesh_p1.msh +gmsh geometry.geo -2 -setnumber element_order 2 -o mesh_p2.msh diff --git a/examples/cooks_membrane/geometry.geo b/examples/cooks_membrane/geometry.geo new file mode 100644 index 0000000..4a7e145 --- /dev/null +++ b/examples/cooks_membrane/geometry.geo @@ -0,0 +1,53 @@ +//////////////////////////////////////////////////////////// +// Cook's membrane (Q1 quadrilateral mesh) +// +// Geometry: +// (0,44) ----------- (48,60) +// | | +// | | +// (0,0) ------------ (48,44) +//////////////////////////////////////////////////////////// + +SetFactory("OpenCASCADE"); + +// Default if nothing is supplied on the command line +DefineConstant[ + element_order = 1 +]; + +// Geometry +Point(1) = {0, 0, 0}; +Point(2) = {48,44, 0}; +Point(3) = {48,60, 0}; +Point(4) = {0,44, 0}; + +Line(1) = {1,2}; +Line(2) = {2,3}; +Line(3) = {3,4}; +Line(4) = {4,1}; + +Curve Loop(1) = {1,2,3,4}; +Plane Surface(1) = {1}; + +// Structured mesh +nx = 32; +ny = 32; + +Transfinite Curve{1,3} = nx + 1; +Transfinite Curve{2,4} = ny + 1; + +Transfinite Surface{1}; +Recombine Surface{1}; + +// Physical groups +Physical Surface("Domain") = {1}; + +Physical Curve("Left") = {4}; +Physical Curve("Right") = {2}; +Physical Curve("Bottom") = {1}; +Physical Curve("Top") = {3}; + +Mesh.ElementOrder = element_order; +Mesh.SecondOrderIncomplete = 0; + +Mesh 2; diff --git a/examples/cooks_membrane/script.jl b/examples/cooks_membrane/script.jl new file mode 100644 index 0000000..dc70f67 --- /dev/null +++ b/examples/cooks_membrane/script.jl @@ -0,0 +1,224 @@ +import FiniteElementContainers as FEC +using FiniteElementContainers +using Gmsh +using StaticArrays +using Tensors + +struct TwoFieldSolidMechanics{NF, NP, NS} <: AbstractPhysics{NF, NP, NS} +end + +struct Displ <: AbstractPhysics{2, 3, 0} +end + +struct Pressure <: AbstractPhysics{1, 3, 0} +end + +function FiniteElementContainers.create_properties(::TwoFieldSolidMechanics) + ρ = 1e3 + K = 1.e9 + G = 1.e6 + return SVector{3, Float64}(ρ, K, G) +end + +function jacobian(∇u) + return det(∇u + one(∇u)) +end + +function pk1_stress_iso(props, ∇u, p) + κ, μ = props[2], props[3] + F = ∇u + one(∇u) + J = det(F) + J_m_13 = 1. / cbrt(J) + J_m_23 = J_m_13 * J_m_13 + I_1 = tr(tdot(F)) + F_inv_T = inv(F)' + P_iso = μ * J_m_23 * (F - (1. / 3.) * I_1 * F_inv_T) + return P_iso +end + +function pk1_stress_vol(props, ∇u, p) + κ, μ = props[2], props[3] + F = ∇u + one(∇u) + J = det(F) + F_inv_T = inv(F)' + # P_vol = 0.5 * κ * (J * J - 1.) * F_inv_T + P_vol = p * J * F_inv_T + return P_vol +end + +material_tangent_iso(props, ∇u, p) = Tensors.gradient(z -> pk1_stress_iso(props, z, p), ∇u) +material_tangent_vol(props, ∇u, p) = Tensors.gradient(z -> pk1_stress_vol(props, z, p), ∇u) + +@inline function FiniteElementContainers.residual( + physics::TwoFieldSolidMechanics, interps, x_el, t, dt, + u_el, u_el_old, state_old_q, state_new_q, props_el +) + u_el, p_el = u_el + interps_u, interps_p = interps + x_el_u, x_el_p = x_el + interps_u = map_interpolants(interps_u, x_el_u) + interps_p = map_interpolants(interps_p, x_el_p) + JxW_u = interps_u.JxW + JxW_p = interps_p.JxW + ∇u_q = interpolate_field_gradients(Displ(), interps_u, u_el) + ∇u_q = modify_field_gradients(PlaneStrain(), ∇u_q) + p_q = interpolate_field_values(Pressure(), interps_p, p_el) + + # constitutive + P_iso = pk1_stress_iso(props_el, ∇u_q, p_q[1]) + P_vol = pk1_stress_vol(props_el, ∇u_q, p_q[1]) + J = jacobian(∇u_q) + + P_q = extract_stress(PlaneStrain(), P_iso + P_vol) + G_q = discrete_gradient(PlaneStrain(), interps_u.∇N_X) + R_u = JxW_u * G_q * P_q + R_p = JxW_p * (J - one(J)) * interps_p.N + + return R_u, R_p +end + +@inline function FiniteElementContainers.stiffness( + physics::TwoFieldSolidMechanics, interps, x_el, t, dt, + u_el, u_el_old, state_old_q, state_new_q, props_el +) + u_el, p_el = u_el + interps_u, interps_p = interps + x_el_u, x_el_p = x_el + interps_u = map_interpolants(interps_u, x_el_u) + interps_p = map_interpolants(interps_p, x_el_p) + JxW_u = interps_u.JxW + JxW_p = interps_p.JxW + ∇u_q = interpolate_field_gradients(Displ(), interps_u, u_el) + ∇u_q = modify_field_gradients(PlaneStrain(), ∇u_q) + p_q = interpolate_field_values(Pressure(), interps_p, p_el) + J_q = jacobian(∇u_q) + F_q = ∇u_q + one(∇u_q) + F_inv_T_q = inv(F_q)' + dPdp_q = extract_stress(PlaneStrain(), J_q * F_inv_T_q) + A_iso = material_tangent_iso(props_el, ∇u_q, p_q[1]) + A_vol = material_tangent_vol(props_el, ∇u_q, p_q[1]) + G_q = discrete_gradient(PlaneStrain(), interps_u.∇N_X) + G_pu_x = J_q .* ( + F_inv_T_q[1, 1] .* interps_u.∇N_X[:, 1] + + F_inv_T_q[1, 2] .* interps_u.∇N_X[:, 2] + ) + G_pu_y = J_q .* ( + F_inv_T_q[2, 1] .* interps_u.∇N_X[:, 1] + + F_inv_T_q[2, 2] .* interps_u.∇N_X[:, 2] + ) + Nd = length(G_pu_x) + + tup = MVector{2 * Nd, eltype(G_pu_x)}(undef) + + for i in 1:Nd + tup[2i-1] = G_pu_x[i] + tup[2i] = G_pu_y[i] + end + + G_pu_q = SVector{2 * Nd, eltype(G_pu_x)}(tup) + K_uu = JxW_u * G_q * extract_stiffness(PlaneStrain(), A_iso + A_vol) * G_q' + K_up = JxW_u * G_q * dPdp_q * interps_p.N' + K_pu = JxW_p * interps_p.N * G_pu_q' + K_pp = zero(SMatrix{length(p_el), length(p_el), Float64, length(p_el)^2}) + return ( + (K_uu, K_up), + (K_pu, K_pp) + ) +end + +mesh_u = UnstructuredMesh(Base.source_dir() * "/mesh_p2.msh") +mesh_p = UnstructuredMesh(Base.source_dir() * "/mesh_p1.msh") + +V_u = FunctionSpace(mesh_u, H1Field, Lagrange) +V_p = FunctionSpace(mesh_p, H1Field, Lagrange) +# V_J = FunctionSpace(mesh_p, H1Field, Lagrange) + +u = VectorFunction(V_u, "displ") +p = ScalarFunction(V_p, "pressure") +# J = ScalarFunction(V_J, "jacobian") + +zero_func(_, _) = 0.0 +displ_func(_, t) = 0.01 * t +dbcs_u = DirichletBC[ + DirichletBC("displ_x", zero_func; nodeset_name = "Left") + DirichletBC("displ_y", zero_func; nodeset_name = "Left") + DirichletBC("displ_x", zero_func; nodeset_name = "Right") + DirichletBC("displ_y", displ_func; nodeset_name = "Right") +] +physics = TwoFieldSolidMechanics{3, 0, 0}() +props = create_properties(physics) +times = TimeStepper(0.0, 1.0, 10) + +dof_u, dof_p = DofManager(u), DofManager(p) +dof = (dof_u, dof_p) +# dof_u, dof_p, dof_J = DofManager(u), DofManager(p), DofManager(J) +# dof = (dof_u, dof_p, dof_J) +asm = FEC.BlockSparseMatrixAssembler(dof) + +p_u = create_parameters(mesh_u, SparseMatrixAssembler(dof_u), physics, props; dirichlet_bcs = dbcs_u, times = times) +p_p = create_parameters(mesh_p, SparseMatrixAssembler(dof_p), physics, props; times = times) +params = (p_u, p_p) + +FEC.update_dofs!( + asm, + (p_u.dirichlet_bcs, p_p.dirichlet_bcs), + (p_u.periodic_bcs, p_p.periodic_bcs) +) + +pp_u = PostProcessor(mesh_u, "u.exo", u) +pp_p = PostProcessor(mesh_p, "p.exo", p) +# solver = NewtonSolver(DirectLinearSolver(asm)) +# integrator = QuasiStaticIntegrator(solver) + +# for n in 1:20 +# evolve!(integrator, params) +# end + +Uu = create_unknowns(asm) +# # U = create_field(asm) + +# assemble_stiffness!(asm, stiffness, Uu, params) +# K = stiffness(asm) + +# K_up = K.blocks[1, 2] +# K_pu = K.blocks[2, 1] + +# temp = K_up .- K_pu' +# display(K_up .- K_pu') + +for n in 1:10 + FiniteElementContainers.update_time!(params) + FiniteElementContainers.update_bc_values!(params, asm) + r0 = -1e6 + for iter in 1:10 + assemble_vector!(asm, residual, Uu, params) + R = residual(asm) + + rnorm = norm(R) + if iter == 1 + r0 = rnorm + end + + if rnorm / r0 < 1e-8 + break + end + + assemble_stiffness!(asm, stiffness, Uu, params) + K = stiffness(asm) + + ΔU = K \ R + Uu .-= ΔU + + println("iter = $iter, |R| = $(rnorm / r0), |ΔU| = $(norm(ΔU))") + + end + + write_times(pp_u, n + 1, params[1].times.time_current) + write_times(pp_p, n + 1, params[2].times.time_current) + write_field(pp_u, n + 1, ("displ_x", "displ_y"), params[1].field) + write_field(pp_p, n + 1, ("pressure",), params[2].field) +end +close(pp_u) +close(pp_p) + +# Δu = K \ R diff --git a/src/Parameters.jl b/src/Parameters.jl index ab64b71..06bc075 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -427,6 +427,21 @@ function update_bc_values!(p::AbstractParameters, assembler) return nothing end +function update_bc_values!(params::Tuple, assembler) + for n in 1:length(params) + p = params[n] + X = coordinates(p) + t = current_time(p) + update_bc_values!(p.dirichlet_bcs, X, t) + # update_bc_values!(p.neumann_bcs, assembler, X, t) + # update_bc_values!(p.periodic_bcs, X, t) + # # update_bc_values!(p.robin_bcs, assembler, X, t, p.field) + # update_source_values!(p.sources, assembler, X, t) + + # update_bc_values!(p[n], assembler) + end +end + function update_bc_values!(p::TypeStableParameters, assembler) X = coordinates(p) t = current_time(p) @@ -454,6 +469,18 @@ function update_dofs!(asm::AbstractAssembler, p::Parameters) return nothing end +function _update_field!(p::AbstractParameters) + p.field_old.data .= p.field.data + return nothing +end + +function _update_field!(p::Tuple) + for n in 1:length(p) + _update_field!(p[n]) + end + return nothing +end + function _update_for_assembly!(p::AbstractParameters, dof::DofManager, Uu) update_field_dirichlet_bcs!(p.field, p.dirichlet_bcs) update_field_unknowns!(p.field, dof, Uu) @@ -477,6 +504,12 @@ function _update_for_assembly!(p::AbstractParameters, dof::DofManager, Uu, Vu) return nothing end +function _update_for_assembly!(p, dof::Tuple, Uu) + for n in 1:length(p) + _update_for_assembly!(p[n], dof[n], view(Uu, BlockArrays.Block(n))) + end +end + # Full-DOF flavor: caller is responsible for assembling the merged # vectors U_full = [Uu; U_BC] and v_full = [v_free; v_BC] themselves. # Unlike the free-DOF flavors above, we do NOT call @@ -501,3 +534,10 @@ function update_time!(p::AbstractParameters) p.times.time_current = current_time(p.times) + time_step(p.times) return nothing end + +function update_time!(p::Tuple) + for n in 1:length(p) + update_time!(p[n]) + end + return nothing +end diff --git a/src/Solvers.jl b/src/Solvers.jl index 608e410..2448cea 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -49,6 +49,17 @@ struct DirectLinearSolver{ # what's the best way to do this with general solvers? ΔUu::U + function DirectLinearSolver(assembler::BlockSparseMatrixAssembler) + preconditioner = I + ΔUu = similar(assembler.residual_unknowns) + fill!(ΔUu, zero(eltype(ΔUu))) + new{typeof(assembler), typeof(preconditioner), typeof(ΔUu)}( + assembler, preconditioner, + DirectLinearSolverSettings(), + TimerOutput(), ΔUu + ) + end + function DirectLinearSolver(assembler::SparseMatrixAssembler) preconditioner = I ΔUu = similar(assembler.residual_unknowns) diff --git a/src/assemblers/Assemblers.jl b/src/assemblers/Assemblers.jl index 9a34da8..e80b209 100644 --- a/src/assemblers/Assemblers.jl +++ b/src/assemblers/Assemblers.jl @@ -2,7 +2,7 @@ $(TYPEDEF) $(TYPEDFIELDS) """ -abstract type AbstractAssembler{Dof <: DofManager} end +abstract type AbstractAssembler end """ $(TYPEDSIGNATURES) """ @@ -79,7 +79,6 @@ function _assemble_element!( for n in axes(conns, 1) global_id = n_dofs * (conns[n] - 1) + d local_id = n_dofs * (n - 1) + d - # Atomix.@atomic storage.data[global_id] += R_el[local_id] fec_atomic_add!(storage, global_id, R_el[local_id]) end end @@ -103,9 +102,6 @@ function _assemble_element!( return nothing end -# TODO we'll need a regular matrix implementation -# as well (Can we live with 1?) -# sparse matrix function _assemble_element!( storage, K_el::SMatrix{NDOF1, NDOF2, T, NDOF1xNDOF2}, conns, # all connectivities for this element @@ -151,7 +147,6 @@ end NNPE = ReferenceFiniteElements.num_cell_dofs(ref_fe) NxNDof = NNPE * NF u_el = @views SMatrix{NF, NNPE, eltype(U), NxNDof}(U[:, conns]) - # u_el = @views SMatrix{NNPE, ND, eltype(U), NxNDof}(U[:, conns]) return u_el end @@ -228,13 +223,17 @@ end return zeros(SVector{NxNDof, eltype(U)}) end -# """ -# $(TYPEDSIGNATURES) -# """ -# function _quadrature_level_state(state::AbstractArray{<:Number, 3}, q::Int, e::Int) -# state_q = view(state, :, q, e) -# return state_q -# end +@inline function _element_scratch( + ::AssembledMatrix, + ref_fe_row, + U_row::H1Field{T, D, NFr}, + ref_fe_col, + U_col::H1Field{T, D, NFc}, +) where {T, D, NFr, NFc} + nr = num_cell_dofs(ref_fe_row) * NFr + nc = num_cell_dofs(ref_fe_col) * NFc + return zeros(SMatrix{nr, nc, T, nr * nc}) +end function _sparse_matrix_mass(asm::AbstractAssembler, coo_storage) type = _sparse_matrix_type(asm) @@ -478,6 +477,7 @@ end include("SparsityPatterns.jl") # types +include("BlockMatrixAssemblers.jl") include("MatrixFreeAssembler.jl") include("SparseMatrixAssembler.jl") diff --git a/src/assemblers/BlockMatrixAssemblers.jl b/src/assemblers/BlockMatrixAssemblers.jl new file mode 100644 index 0000000..6bf1257 --- /dev/null +++ b/src/assemblers/BlockMatrixAssemblers.jl @@ -0,0 +1,261 @@ +abstract type AbstractBlockAssembler <: AbstractAssembler end + +function create_field(asm::AbstractBlockAssembler) + return create_field(asm.dof) +end + +function create_unknowns(asm::AbstractBlockAssembler) + return create_unknowns(asm.dof) +end + +mutable struct BlockSparseMatrixAssembler{ + I <: AbstractVector{Int}, + R <: AbstractVector{Float64}, + D, + F, + U <: BlockedVector, + S +} <: AbstractBlockAssembler + dof::D + matrix_patterns::Matrix{SparseMatrixPattern{I, R}} + vector_patterns::Vector{SparseVectorPattern{I}} + residual_storage::F + residual_unknowns::U + stiffness_storage::S +end + +function BlockSparseMatrixAssembler(dof::Tuple) + matrix_patterns = Matrix{SparseMatrixPattern{Vector{Int}, Vector{Float64}}}(undef, length(dof), length(dof)) + vector_patterns = Vector{SparseVectorPattern{Vector{Int}}}(undef, length(dof)) + for i in 1:length(dof) + for j in 1:length(dof) + matrix_patterns[i, j] = SparseMatrixPattern(dof[i], dof[j]) + end + vector_patterns[i] = SparseVectorPattern(dof[i]) + end + n_matrix_entries = map(num_entries, matrix_patterns) + residual = create_field(dof) + residual_unknowns = create_unknowns(dof) + stiffness_storage = map(zeros, n_matrix_entries) + return BlockSparseMatrixAssembler( + dof, matrix_patterns, vector_patterns, + residual, residual_unknowns, + stiffness_storage + ) +end + +function Base.show(io::IO, asm::BlockSparseMatrixAssembler) + sz = size(asm.matrix_patterns) + println(io, "BlockSparseMatrixAssembler:") + println(io, " Block layout = $(sz[1]) x $(sz[2])") + for dof in asm.dof + show(io, dof; pad = " ") + println(io) + end + println(" Matrix sizes:") + for i in axes(asm.matrix_patterns, 1) + string = " " + for j in axes(asm.matrix_patterns, 2) + string = string * "($(length(asm.dof[i].unknown_dofs)), $(length(asm.dof[j].unknown_dofs)))" + if j < size(asm.matrix_patterns, 2) + string = string * ", " + end + end + println(io, string) + end + # println(io, " Block variables = ") +end + +function assemble_stiffness!( + assembler::BlockSparseMatrixAssembler, func::F, Uu, p +) where F <: Function + @assert length(assembler.dof) == 2 "Only two spaces supported currently" + # storage = assembler.residual_storage + storage = assembler.stiffness_storage + map(x -> fill!(x, zero(eltype(x))), storage) + fspace = map(function_space, assembler.dof) + X = map(coordinates, p) + # should we do a check that all times, and time steps are consistent? + t = current_time(p[1]) + Δt = time_step(p[1]) + U = map(x -> x.field, p) + U_old = map(x -> x.field_old, p) + + for sol_id in 1:length(fspace) + _update_for_assembly!(p[sol_id], assembler.dof[sol_id], Uu[BlockArrays.Block(sol_id)]) + end + + return_type = AssembledMatrix() + conns = map(x -> x.elem_conns.data, fspace) + # conns = map(block_conns, fspace) + coffsets = map(x -> x.elem_conns.offsets, fspace) + physics = p[1].physics + props = p[1].properties + for b in 1:num_blocks(fspace[1]) + block_physics = values(physics)[b] + ref_fe = map(x -> block_reference_element(x, b), fspace) + num_q_pts = map(num_cell_quadrature_points, ref_fe) + @assert all(==(num_q_pts[1]), num_q_pts) + num_q_pts = num_q_pts[1] + for e in 1:block_entity_size(fspace[1], b)[2] + conn = map((r, c, co) -> connectivity(r, c, e, co[b]), ref_fe, conns, coffsets) + out = map((r, c, x, u, u_old) -> element_level_fields(r, c, x, u, u_old), ref_fe, conn, X, U, U_old) + x_el = map(x -> x[1], out) + u_el = map(x -> x[2], out) + u_el_old = map(x -> x[3], out) + props_el = properties(props, e, b) + # val_el = map((r, u) -> _element_scratch(return_type, r, u), ref_fe, U) + nfields = length(U) + + val_el = ntuple(i -> ntuple(j->begin + _element_scratch( + return_type, + ref_fe[i], U[i], + ref_fe[j], U[j] + ) + end, nfields), nfields) + for q in 1:num_q_pts + interps = map(r -> _cell_interpolants(r, q), ref_fe) + state_old_q = state_variables(p[1].state_old, q, e, b) + state_new_q = state_variables(p[1].state_new, q, e, b) + val_q = func(block_physics, interps, x_el, t, Δt, u_el, u_el_old, state_old_q, state_new_q, props_el) + # val_el = map((f, vq, ve) -> _accumulate_q_value(return_type, f, vq, ve, q, e), U, val_q, val_el) + val_el = map( + (vq1, ve1) -> + map((f, vq, ve) -> _accumulate_q_value(return_type, f, vq, ve, q, e), U, vq1, ve1), + val_q, val_el + ) + end + # map((f, v, c) -> _assemble_element!(f, v, c, e), U, val_el, conn, e) + for i in 1:nfields + for j in 1:nfields + _assemble_element!(assembler.stiffness_storage[i, j], val_el[i][j], conn[i], e) + end + end + end + end +end + +function assemble_vector!( + assembler::BlockSparseMatrixAssembler, func::F, Uu, p +) where F <: Function + @assert length(assembler.dof) == 2 "Only two spaces supported currently" + storage = assembler.residual_storage + map(x -> fill!(x, zero(eltype(x))), storage) + fspace = map(function_space, assembler.dof) + X = map(coordinates, p) + # should we do a check that all times, and time steps are consistent? + t = current_time(p[1]) + Δt = time_step(p[1]) + U = map(x -> x.field, p) + U_old = map(x -> x.field_old, p) + + for sol_id in 1:length(fspace) + _update_for_assembly!(p[sol_id], assembler.dof[sol_id], Uu[BlockArrays.Block(sol_id)]) + end + + return_type = AssembledVector() + conns = map(x -> x.elem_conns.data, fspace) + coffsets = map(x -> x.elem_conns.offsets, fspace) + physics = p[1].physics + props = p[1].properties + for b in 1:num_blocks(fspace[1]) + block_physics = values(physics)[b] + ref_fe = map(x -> block_reference_element(x, b), fspace) + num_q_pts = map(num_cell_quadrature_points, ref_fe) + @assert all(==(num_q_pts[1]), num_q_pts) + num_q_pts = num_q_pts[1] + for e in 1:block_entity_size(fspace[1], b)[2] + conn = map((r, c, co) -> connectivity(r, c, e, co[b]), ref_fe, conns, coffsets) + out = map((r, c, x, u, u_old) -> element_level_fields(r, c, x, u, u_old), ref_fe, conn, X, U, U_old) + x_el = map(x -> x[1], out) + u_el = map(x -> x[2], out) + u_el_old = map(x -> x[3], out) + props_el = properties(props, e, b) + val_el = map((r, u) -> _element_scratch(return_type, r, u), ref_fe, U) + for q in 1:num_q_pts + interps = map(r -> _cell_interpolants(r, q), ref_fe) + state_old_q = state_variables(p[1].state_old, q, e, b) + state_new_q = state_variables(p[1].state_new, q, e, b) + val_q = func(block_physics, interps, x_el, t, Δt, u_el, u_el_old, state_old_q, state_new_q, props_el) + val_el = map((f, vq, ve) -> _accumulate_q_value(return_type, f, vq, ve, q, e), U, val_q, val_el) + # @show val_el + end + # out = map((f, v, c) -> _assemble_element!(f, v, c, e), U, val_el, conn) + for sol_id in 1:length(U) + # @show "Hur" + _assemble_element!(assembler.residual_storage[sol_id], val_el[sol_id], conn[sol_id], e) + end + end + # @show assembler.residual_storage + end +end + +# function assemble_vector_neumann_bc!(assembler, Uu, p) +# @warn "Sources not supported in block solvers yet." +# return nothing +# end + +# function assemble_vector_source!(assembler, Uu, p) +# @warn "Sources not supported in block solvers yet." +# return nothing +# end + +function residual(asm::BlockSparseMatrixAssembler) + # @show asm.residual_storage + for (b, (d, s)) in enumerate(zip(asm.dof, asm.residual_storage)) + # @show s + extract_field_unknowns!(view(asm.residual_unknowns, BlockArrays.Block(b)), d, s) + end + # @show asm.residual_unknowns + return asm.residual_unknowns +end + +function stiffness(asm::BlockSparseMatrixAssembler) + ndofs_each = map(x -> length(x.unknown_dofs), asm.dof) |> collect + ndofs = reduce(+, ndofs_each) + K = BlockArray(spzeros(ndofs, ndofs), ndofs_each, ndofs_each) + for i in axes(asm.stiffness_storage, 1) + for j in axes(asm.stiffness_storage, 1) + pat = asm.matrix_patterns[i, j] + vals = asm.stiffness_storage[i, j][pat.unknown_dofs] + temp = sparse(pat.Is, pat.Js, vals) + # display(temp) + # K[BlockArrays.Block(i, j)] = _sparse_matrix_stiffness() + K[BlockArrays.Block(i, j)] = temp + end + end + return K +end + +# this won't work with condensed right now +function update_dofs!( + assembler::BlockSparseMatrixAssembler, dirichlet_bcs, periodic_bcs +) + ddofs = map(dirichlet_dofs, dirichlet_bcs) + pdofs = map(periodic_dofs, periodic_bcs) + pdofs_side_a = map(x -> x[1], pdofs) + pdofs_side_b = map(x -> x[2], pdofs) + + # update dof managers first + for (n, dof) in enumerate(assembler.dof) + update_dofs!(dof, ddofs[n], pdofs_side_a[n], pdofs_side_b[n]) + end + + # now update sparsity patterns + for i in axes(assembler.matrix_patterns, 1) + for j in axes(assembler.matrix_patterns, 2) + _update_dofs!( + assembler.matrix_patterns[i, j], + assembler.dof[i], ddofs[i], pdofs_side_b[i], + assembler.dof[j], ddofs[j], pdofs_side_b[j] + ) + end + end + + # update size of residual unknowns + assembler.residual_unknowns = create_unknowns(assembler.dof) + return nothing +end + +_use_inplace_methods(::BlockSparseMatrixAssembler) = false diff --git a/src/assemblers/MatrixFreeAssembler.jl b/src/assemblers/MatrixFreeAssembler.jl index 1834601..ca0e1ae 100644 --- a/src/assemblers/MatrixFreeAssembler.jl +++ b/src/assemblers/MatrixFreeAssembler.jl @@ -5,7 +5,7 @@ struct MatrixFreeAssembler{ RV <: AbstractArray{Float64, 1}, Var <: AbstractFunction, FieldStorage <: AbstractField{Float64, NumArrDims, RV} -} <: AbstractAssembler{DofManager{Condensed, Int, IV, Var}} +} <: AbstractAssembler dof::DofManager{Condensed, Int, IV, Var} vector_pattern::SparseVectorPattern{IV} constraint_storage::RV diff --git a/src/assemblers/SparseMatrixAssembler.jl b/src/assemblers/SparseMatrixAssembler.jl index 4367a80..8ad555d 100644 --- a/src/assemblers/SparseMatrixAssembler.jl +++ b/src/assemblers/SparseMatrixAssembler.jl @@ -13,7 +13,7 @@ struct SparseMatrixAssembler{ RV <: AbstractArray{Float64, 1}, Var <: AbstractFunction, FieldStorage -} <: AbstractAssembler{DofManager{Condensed, Int, IV, Var}} +} <: AbstractAssembler dof::DofManager{Condensed, Int, IV, Var} matrix_pattern::SparseMatrixPattern{IV, RV} vector_pattern::SparseVectorPattern{IV} diff --git a/src/assemblers/SparsityPatterns.jl b/src/assemblers/SparsityPatterns.jl index 20d9757..051a5c5 100644 --- a/src/assemblers/SparsityPatterns.jl +++ b/src/assemblers/SparsityPatterns.jl @@ -90,8 +90,8 @@ function SparseMatrixPattern(dof_1::DofManager, dof_2::DofManager) conn_2 = connectivity(fspace_2, e, b) dof_conn_1 = @views reshape(ids_1[:, conn_1], ND1 * num_entities_per_element(fspace_1, b)) dof_conn_2 = @views reshape(ids_2[:, conn_2], ND2 * num_entities_per_element(fspace_2, b)) - for i in axes(dof_conn_1, 1) - for j in axes(dof_conn_2, 1) + for j in axes(dof_conn_2, 1) + for i in axes(dof_conn_1, 1) Is[n] = dof_conn_1[i] Js[n] = dof_conn_2[j] unknown_dofs[n] = n @@ -206,10 +206,10 @@ function _update_dofs!( conns_2 = connectivity(fspace_2, e, b) dof_conns_1 = @views reshape(ids_1[:, conns_1], ND1 * num_entities_per_element(fspace_1, b)) dof_conns_2 = @views reshape(ids_2[:, conns_2], ND2 * num_entities_per_element(fspace_2, b)) - for i in axes(dof_conns_1, 1) - ri = dof_to_unknown_index(dof_1, dof_conns_1[i]) - for j in axes(dof_conns_2, 1) - rj = dof_to_unknown_index(dof_2, dof_conns_2[j]) + for j in axes(dof_conns_2, 1) + rj = dof_to_unknown_index(dof_2, dof_conns_2[j]) + for i in axes(dof_conns_1, 1) + ri = dof_to_unknown_index(dof_1, dof_conns_1[i]) if ri > 0 && rj > 0 n_entries += 1 end @@ -235,10 +235,10 @@ function _update_dofs!( conns_2 = connectivity(fspace_2, e, b) dof_conns_1 = @views reshape(ids_1[:, conns_1], ND1 * num_entities_per_element(fspace_1, b)) dof_conns_2 = @views reshape(ids_2[:, conns_2], ND2 * num_entities_per_element(fspace_2, b)) - for i in axes(dof_conns_1, 1) - ri = dof_to_unknown_index(dof_1, dof_conns_1[i]) - for j in axes(dof_conns_2, 1) - rj = dof_to_unknown_index(dof_2, dof_conns_2[j]) + for j in axes(dof_conns_2, 1) + rj = dof_to_unknown_index(dof_2, dof_conns_2[j]) + for i in axes(dof_conns_1, 1) + ri = dof_to_unknown_index(dof_1, dof_conns_1[i]) if ri > 0 && rj > 0 pattern.Is[n] = ri pattern.Js[n] = rj diff --git a/src/integrators/QuasiStaticIntegrator.jl b/src/integrators/QuasiStaticIntegrator.jl index b110b38..7061aef 100644 --- a/src/integrators/QuasiStaticIntegrator.jl +++ b/src/integrators/QuasiStaticIntegrator.jl @@ -21,7 +21,8 @@ function evolve!(integrator::QuasiStaticIntegrator, p) # the call below will ensure the return fields have bcs properly enfroced # before being saved as the old solution for the next step. _update_for_assembly!(p, integrator.solver.linear_solver.assembler.dof, integrator.solution) - p.field_old.data .= p.field.data + # p.field_old.data .= p.field.data + _update_field!(p) # Track convergence: check if the last Newton increment is small. norm_ΔUu = sqrt(sum(abs2, integrator.solver.linear_solver.ΔUu))