diff --git a/src/problems.jl b/src/problems.jl index dd3d275..a303719 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -8,7 +8,7 @@ abstract MixedProblem <: AbstractProblem """ General linearized problem to solve - K*u + C1'*la = f + K*u + C1.T*la = f C2*u + D*la = g """ type Assembly @@ -19,7 +19,7 @@ type Assembly # for boundary assembly C1 :: SparseMatrixCOO C2 :: SparseMatrixCOO - D :: SparseMatrixCOO + D :: SparseMatrixCOO g :: SparseMatrixCOO u :: Vector{Float64} # solution vector u @@ -127,6 +127,24 @@ function initialize!(problem::Problem, time::Real) element[field_name] = (time => data) end end + # if this is boundary problem and not dirichlet problem, initialize field + # for primary variable too + is_boundary_problem(problem) || return + is_dirichlet_problem(problem) && return + field_name = get_parent_field_name(problem) + for element in get_elements(problem) + gdofs = get_gdofs(element, problem) + if haskey(element, field_name) + # if field is found, copy last known solution to new time as initial guess + if !isapprox(last(element[field_name]).time, time) + last_data = copy(last(element[field_name]).data) + push!(element[field_name], time => last_data) + end + else # if field not found at all, initialize new zero field. + data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)] + element[field_name] = (time => data) + end + end end """ Update problem solution vector for assembly. """ @@ -169,20 +187,31 @@ This assumes that element is properly initialized so that last known field data is from current time. For boundary problems solution is updated from lambda vector and for field problems from actual solution vector. """ -function update_elements!(problem, u, la) +function update_elements!{P<:FieldProblem}(problem::Problem{P}, u, la) field_name = get_unknown_field_name(problem) field_dim = get_unknown_field_dimension(problem) nnodes = round(Int, length(u)/field_dim) - - solution = nothing - if is_field_problem(problem) - solution = reshape(u, field_dim, nnodes) - elseif is_boundary_problem(problem) - solution = reshape(la, field_dim, nnodes) - else - error("update_elements!(): unknown problem type $(typeof(problem))") + solution = reshape(u, field_dim, nnodes) + for element in get_elements(problem) + connectivity = get_connectivity(element) # node ids + local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] + last(element[field_name]).data = local_sol end - +end +function update_elements!{P<:BoundaryProblem}(problem::Problem{P}, u, la) + field_name = get_unknown_field_name(problem) + field_dim = get_unknown_field_dimension(problem) + nnodes = round(Int, length(u)/field_dim) + solution = reshape(la, field_dim, nnodes) + for element in get_elements(problem) + connectivity = get_connectivity(element) # node ids + local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] + last(element[field_name]).data = local_sol + end + # if boundary problem is not dirichlet, update also data of main problem + is_dirichlet_problem(problem) && return + field_name = get_parent_field_name(problem) + solution = reshape(u, field_dim, nnodes) for element in get_elements(problem) connectivity = get_connectivity(element) # node ids local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity] @@ -223,4 +252,3 @@ end function push!(problem::Problem, element) push!(problem.elements, element) end - diff --git a/src/solvers.jl b/src/solvers.jl index 58b9de8..80fcb9b 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -362,4 +362,3 @@ function call(solver::Solver) # 3. did not converge throw(NonlinearConvergenceError(solver)) end - diff --git a/src/utils.jl b/src/utils.jl index c4348d3..0825979 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -15,16 +15,18 @@ field_dim degrees of freedom / node elements elements used to calculate vector +vec_dim + used to resize solution vector if given time """ -function calculate_nodal_vector(field_name::ASCIIString, field_dim::Int, - elements::Vector{Element}, time::Real) +function calculate_nodal_vector(field_name, field_dim, elements::Vector{Element}, + time, vec_dim=0) A = SparseMatrixCOO() b = SparseMatrixCOO() for element in elements haskey(element, field_name) || continue gdofs = get_gdofs(element, 1) - for ip in get_integration_points(element, Val{2}) + for ip in get_integration_points(element, Val{3}) J = get_jacobian(element, ip, time) w = ip.weight*norm(J) f = element(field_name, ip, time) @@ -40,11 +42,18 @@ function calculate_nodal_vector(field_name::ASCIIString, field_dim::Int, nz = sort(unique(rowvals(A))) x = zeros(size(b)...) x[nz, :] = A[nz,nz] \ b[nz, :] - return vec(transpose(x)) + x = vec(transpose(x)) + if vec_dim != 0 + v = zeros(vec_dim) + v[1:length(x)] = x + return v + else + return x + end end -function calculate_rotated_nodal_vector(field_name::ASCIIString, field_dim::Int, - elements::Vector{Element}, time::Real) +function calculate_rotated_nodal_vector(field_name, field_dim, elements::Vector{Element}, + time, vec_dim=0) A = SparseMatrixCOO() b = SparseMatrixCOO() for element in elements @@ -68,7 +77,13 @@ function calculate_rotated_nodal_vector(field_name::ASCIIString, field_dim::Int, nz = sort(unique(rowvals(A))) x = zeros(size(b)...) x[nz, :] = A[nz,nz] \ b[nz, :] - return vec(transpose(x)) + if vec_dim != 0 + v = zeros(vec_dim) + v[1:length(x)] = x + return v + else + return x + end end """ Collect normal-tangential coordinates to rotation matrix Q. @@ -93,4 +108,3 @@ function get_rotation_matrix(elements, time) end return R end -