diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 7d679d3..f33bd2c 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -60,7 +60,7 @@ include("assembly.jl") include("solver_utils.jl") include("solvers.jl") export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver, - get_unknown_field_name, get_formulation_type, + get_unknown_field_name, get_formulation_type, get_problems, get_field_problems, get_boundary_problems, get_field_assembly, get_boundary_assembly, initialize!, create_projection, eliminate_interior_dofs @@ -86,6 +86,7 @@ export calculate_normals, include("problems_contact.jl") include("problems_contact_2d.jl") include("problems_contact_3d.jl") +include("problems_contact_2d_autodiff.jl") export Contact module API diff --git a/src/problems.jl b/src/problems.jl index 4c76620..ae5cf18 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -203,7 +203,7 @@ function update_assembly!(problem, u, la) end # copy current solutions to previous ones and add/replace new solution - # TODO: here we have couple of options and they needs to be clarified + # TODO: here we have couple of options and they need to be clarified # for total formulation we are solving total quantity Ku = f while in # incremental formulation we solve KΔu = f and u = u + Δu assembly.u_prev = copy(assembly.u) @@ -217,7 +217,7 @@ function update_assembly!(problem, u, la) assembly.u += u assembly.la = la elseif get_formulation_type(problem) == :forwarddiff - info("$(problem.name): forwarddiff formulation, adding increment to solution vector") + info("$(problem.name): forwarddiff formulation, adding increment to solution vector and reaction force vector") assembly.u += u assembly.la += la else diff --git a/src/problems_contact.jl b/src/problems_contact.jl index 08d0fec..0316fbc 100644 --- a/src/problems_contact.jl +++ b/src/problems_contact.jl @@ -34,7 +34,11 @@ function get_unknown_field_name(problem::Problem{Contact}) end function get_formulation_type(problem::Problem{Contact}) - return :incremental + if problem.properties.use_forwarddiff + return :forwarddiff + else + return :incremental + end end function assemble!(problem::Problem{Contact}, time::Real) @@ -50,3 +54,4 @@ function assemble!(problem::Problem{Contact}, time::Real) assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff) end +typealias ContactElements2D Union{Seg2} diff --git a/src/problems_contact_2d.jl b/src/problems_contact_2d.jl index 52b505a..20fc8a0 100644 --- a/src/problems_contact_2d.jl +++ b/src/problems_contact_2d.jl @@ -1,17 +1,10 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -typealias ContactElements2D Union{Seg2} - """ Frictionless 2d small sliding contact without forwarddiff. -problem -time -dimension -finite_sliding -friction -use_forwarddiff +true/false flags: finite_sliding, friction, use_forwarddiff """ function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}}, diff --git a/src/problems_contact_2d_autodiff.jl b/src/problems_contact_2d_autodiff.jl index 7849b73..68b1608 100644 --- a/src/problems_contact_2d_autodiff.jl +++ b/src/problems_contact_2d_autodiff.jl @@ -1,6 +1,8 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +using ForwardDiff + """ Find segment from slave element corresponding to master element nodes. Parameters @@ -21,10 +23,10 @@ function project_from_master_to_slave{E<:MortarElements2D}( slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector; tol=1.0e-10, max_iterations=20) - x1(xi1) = vec(get_basis(E, xi1))*x1_ - dx1(xi1) = vec(get_dbasis(E, xi1))*x1_ - n1(xi1) = vec(get_basis(E, xi1))*n1_ - dn1(xi1) = vec(get_dbasis(E, xi1))*n1_ + x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_ + dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_ + n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_ + dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_ cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) @@ -53,8 +55,8 @@ function project_from_slave_to_master{E<:MortarElements2D}( master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI; tol=1.0e-10, max_iterations=20) - x2(xi2) = vec(get_basis(E, xi2))*x2_ - dx2(xi2) = vec(get_dbasis(E, xi2))*x2_ + x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_ + dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_ cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi2) = cross2(x2(xi2)-x1, n1) dR(xi2) = cross2(dx2(xi2), n1) @@ -73,12 +75,19 @@ function project_from_slave_to_master{E<:MortarElements2D}( end -""" Assemble Mortar problem for two-dimensional problems, i.e. for Seg2 and Seg3 elements. """ -function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) +""" +Frictionless 2d finite sliding contact with forwarddiff. + +true/false flags: finite_sliding, friction, use_forwarddiff +""" +function assemble!(problem::Problem{Contact}, time::Float64, + ::Type{Val{1}}, ::Type{Val{true}}, + ::Type{Val{false}}, ::Type{Val{true}}) props = problem.properties field_dim = get_unknown_field_dimension(problem) field_name = get_parent_field_name(problem) + slave_elements = get_slave_elements(problem) function calculate_interface(x::Vector) @@ -94,16 +103,15 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) # 1. update nodal normals for slave elements Q = [0.0 -1.0; 1.0 0.0] normals = zeros(u) - for element in get_elements(problem) - haskey(element, "master elements") || continue + for element in slave_elements conn = get_connectivity(element) push!(S, conn...) gdofs = get_gdofs(element, field_dim) X_el = element("geometry", time) u_el = Field(Vector[u[:,i] for i in conn]) x_el = X_el + u_el - for ip in get_integration_points(element, Val{3}) - dN = get_dbasis(element, ip) + for ip in get_integration_points(element, 3) + dN = get_dbasis(element, ip, time) N = element(ip, time) t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]) normals[:, conn] += ip.weight*Q*t'*N @@ -118,10 +126,14 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) normals[:,i] = -normals[:,i] end end + normals2 = Dict() + for j in S + normals2[j] = normals[:,j] + end + update!(slave_elements, "normal", time => normals2) # 2. loop all slave elements - for slave_element in get_elements(problem) - haskey(slave_element, "master elements") || continue + for slave_element in slave_elements slave_element_nodes = get_connectivity(slave_element) X1 = slave_element("geometry", time) @@ -130,20 +142,19 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) la1 = Field(Vector[la[:,i] for i in slave_element_nodes]) n1 = Field(Vector[normals[:,i] for i in slave_element_nodes]) nnodes = size(slave_element, 2) - update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data)) # 3. loop all master elements - for master_element in slave_element["master elements"] + for master_element in slave_element("master elements", time) master_element_nodes = get_connectivity(master_element) X2 = master_element("geometry", time) u2 = Field(Vector[u[:,i] for i in master_element_nodes]) x2 = X2 + u2 - x1_midpoint = 1/2*(x1[1]+x1[2]) - x2_midpoint = 1/2*(x2[1]+x2[2]) - distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint)) - distance > props.maximum_distance && continue + #x1_midpoint = 1/2*(x1[1]+x1[2]) + #x2_midpoint = 1/2*(x2[1]+x2[2]) + #distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint)) + #distance > props.maximum_distance && continue # calculate segmentation: we care only about endpoints # note: these are quadratic/cubic functions, analytical solution possible @@ -151,7 +162,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) xi1b = -Inf try xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1]) - xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end]) + xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2]) catch info("failed to create projection!!!!") # TODO @@ -163,13 +174,14 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) De = zeros(nnodes, nnodes) Me = zeros(nnodes, nnodes) - for ip in get_integration_points(slave_element, Val{5}) + for ip in get_integration_points(slave_element, 3) # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip) + dN = get_dbasis(slave_element, ip, time) j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) w = ip.weight*norm(j)*l - xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N1 = get_basis(slave_element, xi_s) + xi = ip.coords[1] + xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) + N1 = get_basis(slave_element, xi_s, time) De += w*diagm(vec(N1)) Me += w*N1'*N1 end @@ -179,25 +191,26 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) master_dofs = get_gdofs(master_element, field_dim) # 4. loop integration points of segment - for ip in get_integration_points(slave_element, Val{5}) + for ip in get_integration_points(slave_element, 3) # jacobian of slave element in deformed state - dN = get_dbasis(slave_element, ip) + dN = get_dbasis(slave_element, ip, time) j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)]) w = ip.weight*norm(j)*l # project gauss point from slave element to master element - xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1) - N1 = vec(get_basis(slave_element, xi_s)) + xi = ip.coords[1] + xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1) + N1 = vec(get_basis(slave_element, xi_s, time)) x_s = N1*x1 # coordinate in gauss point n_s = N1*n1 # normal direction in gauss point t_s = Q'*n_s # tangent direction in gauss point xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2) - N2 = vec(get_basis(master_element, xi_m)) + N2 = vec(get_basis(master_element, xi_m, time)) x_m = N2*x2 Phi = Ae*N1 la_s = Phi*la1 # traction force in gauss point - gn = props.gap_sign*dot(n_s, x_s - x_m) # normal gap + gn = -dot(n_s, x_s - x_m) # normal gap fc[:,slave_element_nodes] += w*la_s*N1' fc[:,master_element_nodes] -= w*la_s*N2' @@ -213,22 +226,22 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) # at this point we have calculated contact force fc and gap for all slave elements. # next task is to find out are they in contact or not and remove inactive nodes - nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap)))) - info("gap: $nzgap") + #nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap)))) + #info("gap: $nzgap") for (i, j) in enumerate(sort(collect(S))) - if j in props.always_inactive - info("special node $j always inactive") - C[:,j] = la[:,j] - continue - end +# if j in props.always_inactive +# info("special node $j always inactive") +# C[:,j] = la[:,j] +# continue +# end n = normals[:,j] t = Q'*n lan = dot(n, la[:,j]) lat = dot(t, la[:,j]) if lan - gap[1, j] > 0 - info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))") +# info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))") C[1,j] += gap[1, j] C[2,j] += lat else @@ -242,22 +255,23 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) # x doesn't mean deformed configuration here x = [problem.assembly.u; problem.assembly.la] - ndofs = round(Int, length(x)/2) - A, allresults = ForwardDiff.jacobian(calculate_interface, x, - ForwardDiff.AllResults, cache=autodiffcache) - b = -ForwardDiff.value(allresults) - + if length(x) == 0 + error("2d autodiff contact problem: initialize problem.assembly.u & la before solution") + end + A = ForwardDiff.jacobian(calculate_interface, x) + b = calculate_interface(x) A = sparse(A) b = sparse(b) SparseMatrix.droptol!(A, 1.0e-12) SparseMatrix.droptol!(b, 1.0e-12) + ndofs = round(Int, length(x)/2) K = A[1:ndofs,1:ndofs] C1 = transpose(A[1:ndofs,ndofs+1:end]) C2 = A[ndofs+1:end,1:ndofs] D = A[ndofs+1:end,ndofs+1:end] - f = b[1:ndofs] - g = b[ndofs+1:end] + f = -b[1:ndofs] + g = -b[ndofs+1:end] empty!(problem.assembly) add!(problem.assembly.K, K) @@ -267,6 +281,5 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}) add!(problem.assembly.f, f) add!(problem.assembly.g, g) - return problem.assembly - end + diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl index d4d2284..4e859d1 100644 --- a/src/problems_mortar.jl +++ b/src/problems_mortar.jl @@ -21,14 +21,11 @@ function get_unknown_field_name(problem::Problem{Mortar}) end function get_formulation_type(problem::Problem{Mortar}) - return :incremental - #= if problem.properties.use_forwarddiff return :forwarddiff else return :incremental end - =# end function assemble!(problem::Problem{Mortar}, time::Float64) diff --git a/src/solvers.jl b/src/solvers.jl index c69bf09..3906733 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -337,9 +337,33 @@ end """ Default initializer for solver. """ function initialize!(solver::Solver; show_info=true) show_info && info("Initializing problems ...") + problems = get_problems(solver) + length(problems) != 0 || error("Empty solver, add problems to solver using push!") t0 = Base.time() - for problem in solver.problems + field_problems = get_field_problems(solver) + length(field_problems) != 0 || warn("No field problem found from solver, add some..?") + field_dim = get_unknown_field_dimension(first(field_problems)) + field_name = get_unknown_field_name(first(field_problems)) + info("initialize!(): looks we are solving $field_name, $field_dim dofs/node") + nodes = Set{Int64}() + for problem in problems initialize!(problem, solver.time) + for element in get_elements(problem) + conn = get_connectivity(element) + push!(nodes, conn...) + end + end + nnodes = length(nodes) + info("Total number of nodes in problems: $nnodes") + maxdof = maximum(nnodes)*field_dim + info("# of max dof (=size of solution vector) is $maxdof") + u = zeros(maxdof) + la = zeros(maxdof) + # TODO: this could be used to initialize elements too... + for problem in problems + problem.assembly.u = u + problem.assembly.la = la + # initialize(problem, ....) end t1 = round(Base.time()-t0, 2) show_info && info("Initialized problems in $t1 seconds.") diff --git a/test/test_contact_2d_finite_sliding.jl b/test/test_contact_2d_finite_sliding.jl new file mode 100644 index 0000000..8e35a7f --- /dev/null +++ b/test/test_contact_2d_finite_sliding.jl @@ -0,0 +1,56 @@ +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Postprocess +using JuliaFEM.Test + +@testset "2d curved block with frictionless finite sliding contact using forwarddiff" begin + # FIXME: needs verification of some other fem software + meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d_curved.med" + mesh = aster_read_mesh(meshfile) + + upper = Problem(Elasticity, "upper", 2) + upper.properties.formulation = :plane_stress + upper.properties.finite_strain = true + upper.properties.geometric_stiffness = true + upper.elements = create_elements(mesh, "UPPER") + update!(upper, "youngs modulus", 96.0) + update!(upper, "poissons ratio", 1/3) + + lower = Problem(Elasticity, "lower", 2) + lower.properties.formulation = :plane_stress + lower.properties.finite_strain = true + lower.properties.geometric_stiffness = true + lower.elements = create_elements(mesh, "LOWER") + update!(lower, "youngs modulus", 96.0) + update!(lower, "poissons ratio", 1/3) + + bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement") + bc_upper.elements = create_elements(mesh, "UPPER_TOP") + update!(bc_upper, "displacement 1", 0.0) + update!(bc_upper, "displacement 2", -0.15) + + bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement") + bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") + update!(bc_lower, "displacement 1", 0.0) + update!(bc_lower, "displacement 2", 0.0) + + contact = Problem(Contact, "contact between upper and lower block", 2, "displacement") + contact.properties.rotate_normals = true + contact.properties.finite_sliding = true + contact.properties.friction = false + contact.properties.use_forwarddiff = true + contact_slave_elements = create_elements(mesh, "LOWER_TOP") + contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") + update!(contact_slave_elements, "master elements", contact_master_elements) + contact.elements = [contact_master_elements; contact_slave_elements] + + solver = NonlinearSolver(upper, lower, bc_upper, bc_lower, contact) + solver() + normu = norm(contact.assembly.u) + info("displacement vector norm = $normu") + + # while accurate solution is unknown this is very close to linear solution + # sqrt( ((Stress 11 - Stress 22)^2 + (Stress 22 - Stress 33)^2 + (Stress 33-Stress 11)^2 + 6*(Stress 12^2 + Stress 23^2 + Stress 13^2))/2 ) +# @test isapprox(normu, 0.49745873784105105) + @test isapprox(normu, 0.49745872893844145) +end diff --git a/test/test_contact_2d_small_sliding.jl b/test/test_contact_2d_small_sliding.jl index 9495fcd..5ea82f8 100644 --- a/test/test_contact_2d_small_sliding.jl +++ b/test/test_contact_2d_small_sliding.jl @@ -38,16 +38,16 @@ function get_model(::Type{Val{Symbol("curved 2d contact small sliding")}}) update!(bc_lower, "displacement 1", 0.0) update!(bc_lower, "displacement 2", 0.0) - interface = Problem(Contact, "contact between upper and lower block", 2, "displacement") - interface.properties.dimension = 1 - interface.properties.rotate_normals = true - interface_slave_elements = create_elements(mesh, "LOWER_TOP") - interface_master_elements = create_elements(mesh, "UPPER_BOTTOM") - update!(interface_slave_elements, "master elements", interface_master_elements) - interface.elements = [interface_master_elements; interface_slave_elements] + contact = Problem(Contact, "contact between upper and lower block", 2, "displacement") + contact.properties.dimension = 1 + contact.properties.rotate_normals = true + contact_slave_elements = create_elements(mesh, "LOWER_TOP") + contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") + update!(contact_slave_elements, "master elements", contact_master_elements) + contact.elements = [contact_master_elements; contact_slave_elements] solver = Solver(Nonlinear) - push!(solver, upper, lower, bc_upper, bc_lower, interface) + push!(solver, upper, lower, bc_upper, bc_lower, contact) return solver end @@ -56,23 +56,13 @@ end # FIXME: needs verification of some other fem software solver = get_model("curved 2d contact small sliding") solver() - upper, lower, bc_upper, bc_lower, interface = solver.problems - @test isapprox(norm(interface.assembly.u), 0.49563347601324315) -end - - -function get_mesh(::Type{Val{Symbol("hertz contact, full 2d model")}}) - meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/hertz_2d_full.med" - mesh = aster_read_mesh(meshfile) + upper, lower, bc_upper, bc_lower, contact = solver.problems + @test isapprox(norm(contact.assembly.u), 0.49563347601324315) end function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) - # from fenet d3613 advanced finite element contact benchmarks - # a = 6.21 mm, pmax = 3585 MPa - # this is a very sparse mesh and for that reason pmax is not very - # (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact - # instead integrate pressure in normal and tangential direction - mesh = get_mesh("hertz contact, full 2d model") + meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/hertz_2d_full.med" + mesh = aster_read_mesh(meshfile) upper = Problem(Elasticity, "CYLINDER", 2) upper.properties.formulation = :plane_strain @@ -106,6 +96,9 @@ function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) contact = Problem(Contact, "contact between block and cylinder", 2, "displacement") contact.properties.rotate_normals = true + contact.properties.finite_sliding = false + contact.properties.friction = false + contact.properties.use_forwarddiff = false contact_slave_elements = create_elements(mesh, "CYLINDER_TO_BLOCK") contact_master_elements = create_elements(mesh, "BLOCK_TO_CYLINDER") update!(contact_slave_elements, "master elements", contact_master_elements) @@ -118,14 +111,21 @@ function get_model(::Type{Val{Symbol("hertz contact, full 2d model")}}) end @testset "test frictionless hertz contact, 2d plane strain" begin + # from fenet d3613 advanced finite element contact benchmarks + # a = 6.21 mm, pmax = 3585 MPa + # this is a very sparse mesh and for that reason pmax is not very + # (only 6 elements in -20 .. 20 mm contact zone, 3 elements in contact + # instead integrate pressure in normal and tangential direction solver = get_model("hertz contact, full 2d model") - solver() upper, lower, bc_fixed, bc_sym_23, load, contact = solver.problems + solver() slaves = get_slave_elements(contact) node_ids, la = get_nodal_vector(slaves, "reaction force", 0.0) node_ids, n = get_nodal_vector(slaves, "normal", 0.0) pres = [dot(ni, lai) for (ni, lai) in zip(n, la)] - @test isapprox(maximum(pres), 4060.010799583303) + #@test isapprox(maximum(pres), 4060.010799583303) + # 12 % error in maximum pressure + @test isapprox(maximum(pres), 3585.0; rtol = 12.0e-2) # integrate pressure in normal and tangential direction Rn = 0.0 Rt = 0.0 @@ -141,7 +141,8 @@ end Rt += w*dot(t, la) end end - @test isapprox(Rn, 35.0e3; rtol=0.0015) + # under 0.15 % error in reaction force + @test isapprox(Rn, 35.0e3; rtol=0.15e-2) @test isapprox(Rt, 0.0; atol=10.0) end diff --git a/test/test_heat.jl b/test/test_heat.jl index d3bae0c..db47281 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -39,7 +39,7 @@ using JuliaFEM.Postprocess @test isapprox(T, T_expected; rtol=1.0e-6) end -@testset "one element heat problem" begin +@testset "2d heat problem (one element)" begin X = Dict{Int, Vector{Float64}}( 1 => [0.0,0.0], @@ -65,15 +65,9 @@ end problem.properties.formulation = "2D" push!(problem, el1, el2) - # define boundary element for dirichlet boundary condition - el3 = Element(Seg2, [3, 4]) - update!(el3, "geometry", X) - update!(el3, "temperature 1", 0.0) - - boundary_condition = Problem(Dirichlet, "T=0 on top", 1, "temperature") - push!(boundary_condition, el3) - - # manual assembling of problem + solution: + # Set constant source f=12 with k=6. Accurate solution is + # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x) + # when boundary flux not active (at t=0) assemble!(problem, 0.0) A = full(problem.assembly.K) b = full(problem.assembly.f) @@ -86,26 +80,14 @@ end @test isapprox(A, A_expected) @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0]) - # using Solver - solver = LinearSolver("solve heat problem") - push!(solver, problem, boundary_condition) - - # Set constant source f=12 with k=6. Accurate solution is - # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x) - # when boundary flux not active (at t=0) - solver.time = 0.0 - solver() - # interpolate temperature at middle of element 2 (flux boundary) at time t=0: - T = el2("temperature", [0.0], 0.0) - @test isapprox(T[1], 1.0) - # Set constant flux g=6 on boundary. Accurate solution is # u(x,y) = x which equals T=1 on boundary. # at time t=1.0 all loads should be on. - solver.time = 1.0 - solver() - T = el2("temperature", [0.0], 1.0) - @test isapprox(T[1], 2.0) + empty!(problem) + assemble!(problem, 1.0) + A = full(problem.assembly.K) + b = full(problem.assembly.f) + @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0]) end function T_acc(x)