From ccf32d5b76a76aa88d8f97ed4e87e6711840e9bd Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 13 Dec 2015 13:47:45 +0200 Subject: [PATCH] fixed problems with surface loads. det(element, ip, time) should be avoided. --- src/elasticity.jl | 29 +++++++++++++++-------------- src/elements.jl | 30 ++++++++++++++++++++++++++---- src/equations.jl | 16 +++++++--------- src/heat.jl | 11 ++++++----- src/linear_elasticity.jl | 22 ++++++++++++---------- test/test_elasticity.jl | 12 ++++++------ test/test_elements.jl | 13 ++++++++++++- test/test_linear_elasticity.jl | 18 +++++++++--------- test/test_potential_energy.jl | 9 ++++++--- test/test_solver.jl | 5 +++-- 10 files changed, 102 insertions(+), 63 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index 83db89e..29f59cc 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -70,20 +70,20 @@ https://en.wikipedia.org/wiki/Hooke's_law """ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing) - - # u = element("displacement", ip, time, variation) - r = zeros(Float64, problem.dim, length(element)) + J = get_jacobian(element, ip, time) + # internal forces if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") u = element("displacement", time, variation) grad = element(ip, time, Val{:grad}) gradu = grad*u - F = I + gradu # deformation gradient -# info("gradu = \n$(ForwardDiff.get_value(gradu))") + # deformation gradient + F = I + gradu + # material young = element("youngs modulus", ip, time) poisson = element("poissons ratio", ip, time) mu = young/(2*(1+poisson)) @@ -91,29 +91,30 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element: if P == PlaneStressElasticityProblem lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems end - E = 1/2*(F'*F - I) # strain + + # strain + E = 1/2*(F'*F - I) + + # stress S = lambda*trace(E)*I + 2*mu*E - #J = det(element, ip, time) - #T = J^-1*F*S*F' - #ip["cauchy stress"] = T - #ip["gl strain"] = E - - r += F*S*grad + r += F*S*grad*det(J) end # external forces - volume load if haskey(element, "displacement load") basis = element(ip, time) b = element("displacement load", ip, time) - r -= b*basis + r -= b*basis*det(J) end # external forces - surface traction force if haskey(element, "displacement traction force") basis = element(ip, time) T = element("displacement traction force", ip, time) - r -= T*basis + JT = transpose(J) + s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2]) + r -= T*basis*norm(s) end return vec(r) diff --git a/src/elements.jl b/src/elements.jl index 2ea6884..1510562 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -197,13 +197,35 @@ function call(element::Element, field_name::ASCIIString) return element[field_name] end -function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Number=0.0) + +""" Return the jacobian of element. """ +function get_jacobian{E}(element::Element{E}, xi::Vector{Float64}, time::Real) X = element("geometry", time) - dN = get_dbasis(E, ip.xi) + dN = get_dbasis(E, xi) J = sum([kron(dN[:,i], X[i]') for i=1:length(X)]) - m, n = size(J) - return m == n ? det(J) : norm(J) + return J end +function get_jacobian{E}(element::Element{E}, ip::IntegrationPoint, time::Real) + return get_jacobian(element, ip.xi, time) +end + +""" Return the determinant of jacobian. """ +function LinAlg.det{E<:AbstractElement}(element::Element{E}, xi::Vector{Float64}, time::Real) + warn("det(element, ip, time) is ambiguous: use J = get_jacobian(element, ip, time); det(J) instead.") + J = get_jacobian(element, xi, time) + n, m = size(J) + if n == m + return det(J) + end + JT = transpose(J) + s = size(JT, 2) == 1 ? norm(JT) : norm(cross(JT[:,1], JT[:,2])) + return s +end +function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Real) + return det(element, ip.xi, time) +end + + """ Check does field exist. """ function Base.haskey(element::Element, what) diff --git a/src/equations.jl b/src/equations.jl index 186a538..f284669 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -77,15 +77,15 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time: # 1. if equations are defined we just integrate them, without caring how they are done if has_mass_matrix(problem, element) || has_stiffness_matrix(problem, element) || has_force_vector(problem, element) for ip in get_integration_points(element) - s = ip.weight*det(element, ip, time) + w = ip.weight*det(J) if has_mass_matrix(element) - add!(assembly.mass_matrix, gdofs, gdofs, s*get_mass_matrix(problem, element, ip, time)) + add!(assembly.mass_matrix, gdofs, gdofs, w*get_mass_matrix(problem, element, ip, time)) end if has_stiffness_matrix(element) - add!(assembly.stiffness_matrix, gdofs, gdofs, s*get_stiffness_matrix(problem, element, ip, time)) + add!(assembly.stiffness_matrix, gdofs, gdofs, w*get_stiffness_matrix(problem, element, ip, time)) end if has_force_vector(element) - add!(assembly.force_vector, gdofs, s*get_force_vector(problem, element, ip, time)) + add!(assembly.force_vector, gdofs, w*get_force_vector(problem, element, ip, time)) end end # external loads -- if any nodal loads is defined add to force vector @@ -104,16 +104,15 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time: df = similar(field, data) # integrate potential energy for ip in get_integration_points(element) - s = ip.weight*det(element, ip, time) dw = get_potential_energy(problem, element, ip, time; variation=df) - W += s*dw + W += ip.weight*dw end # external energy -- if any nodal loads is defined, decrease from potential energy if haskey(element, "$unknown_field_name nodal load") P = element["$unknown_field_name nodal load"](time) W -= dot(vec(P), vec(df)) end - return isa(W, Array) ? W[1] : W + return W[1] end hessian, allresults = ForwardDiff.hessian(calc_W, vec(field), AllResults, cache=autodiffcache) @@ -133,9 +132,8 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time: gauss_fields = IntegrationPoint[] # integrate residual vector for ip in get_integration_points(element) - s = ip.weight*det(element, ip, time) dr = get_residual_vector(problem, element, ip, time; variation=df) - R += s*dr + R += ip.weight*dr if ip.changed push!(gauss_fields, ip) end diff --git a/src/heat.jl b/src/heat.jl index 3f1537b..4e34743 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -46,24 +46,25 @@ function assemble!(assembly::Assembly, problem::Problem{HeatProblem}, element::E gdofs = get_gdofs(element, problem.dim) for ip in get_integration_points(element) - w = ip.weight*det(element, ip, time) + w = ip.weight + J = get_jacobian(element, ip, time) N = element(ip, time) if haskey(element, "density") rho = element("density", ip, time) - add!(assembly.mass_matrix, gdofs, gdofs, w*rho*N'*N) + add!(assembly.mass_matrix, gdofs, gdofs, w*rho*N'*N*det(J)) end if haskey(element, "temperature thermal conductivity") dN = element(ip, time, Val{:grad}) k = element("temperature thermal conductivity", ip, time) - add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN) + add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN*det(J)) end if haskey(element, "temperature load") f = element("temperature load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*f) + add!(assembly.force_vector, gdofs, w*N'*f*det(J)) end if haskey(element, "temperature flux") g = element("temperature flux", ip, time) - add!(assembly.force_vector, gdofs, w*N'*g) + add!(assembly.force_vector, gdofs, w*N'*g*norm(J)) end end end diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl index 40d2dc4..629efdc 100644 --- a/src/linear_elasticity.jl +++ b/src/linear_elasticity.jl @@ -16,7 +16,8 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble ndim, nnodes = size(E) B = zeros(6, 3*nnodes) for ip in get_integration_points(element) - w = ip.weight*det(element, ip, time) + w = ip.weight + J = get_jacobian(element, ip, time) N = element(ip, time) if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") v = element("poissons ratio", ip, time) @@ -44,16 +45,16 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble B[6, 3*(i-1)+1] = dN[3,i] B[6, 3*(i-1)+3] = dN[1,i] end - add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B) + add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) end if haskey(element, "displacement load") b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b) + add!(assembly.force_vector, gdofs, w*N'*b*det(J)) end if haskey(element, "displacement traction force") T = element("displacement traction force", ip, time) - L = w*T*N -# dump(L) + JT = transpose(J) + L = w*T*N*norm(cross(JT[:,1], JT[:,2])) add!(assembly.force_vector, gdofs, vec(L)) end end @@ -72,7 +73,8 @@ function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assem ndim, nnodes = size(E) B = zeros(3, 2*nnodes) for ip in get_integration_points(element) - w = ip.weight*det(element, ip, time) + w = ip.weight + J = get_jacobian(element, ip, time) N = element(ip, time) if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") nu = element("poissons ratio", ip, time) @@ -89,17 +91,17 @@ function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assem B[3, 2*(i-1)+1] = dN[2,i] B[3, 2*(i-1)+2] = dN[1,i] end - add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B) + add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) end if haskey(element, "displacement load") b = element("displacement load", ip, time) - add!(assembly.force_vector, gdofs, w*N'*b) + add!(assembly.force_vector, gdofs, w*N'*b*det(J)) end if haskey(element, "displacement traction force") T = element("displacement traction force", ip, time) - L = w*T*N -# dump(L) + L = w*T*N*norm(J) add!(assembly.force_vector, gdofs, vec(L)) end end end + diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 36c9268..7ba22e0 100644 --- a/test/test_elasticity.jl +++ b/test/test_elasticity.jl @@ -92,7 +92,7 @@ function test_continuum_elasticity_with_surface_load() set_geometry!(element1, nodes) # element1["youngs modulus"] = 900.0 # element1["poissons ratio"] = 0.25 - element1["youngs modulus"] = 9000.0 + element1["youngs modulus"] = 900.0 element1["poissons ratio"] = 0.25 element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8]) @@ -105,7 +105,6 @@ function test_continuum_elasticity_with_surface_load() push!(problem, element1) push!(problem, element2) - #= free_dofs = zeros(Bool, 8, 3) x = 1 y = 2 @@ -126,8 +125,8 @@ function test_continuum_elasticity_with_surface_load() info("initial stiffness matrix") dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs]) solve!(problem, free_dofs, 0.0; max_iterations=10) - =# +#= dx = Quad4([1, 4, 8, 5]) dx["displacement 1"] = 0.0 dy = Quad4([1, 5, 6, 2]) @@ -145,17 +144,18 @@ function test_continuum_elasticity_with_surface_load() solver.dump_matrices = true solver.name = "3d_hex8" solver(0.0) +=# disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) info("displacement at tip: $disp") info("displacement on element: ") for (i, d) in enumerate(element1("displacement", 0.0)) - @printf "%d %f %f %f\n" [i;d]... + @printf "%d % f % f % f\n" [i;d]... end # verified using Code Aster. # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm - # @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) - @test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02]) + @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) + #@test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02]) end #test_continuum_elasticity_with_surface_load() diff --git a/test/test_elements.jl b/test/test_elements.jl index 6930d39..a17e3ff 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -6,7 +6,7 @@ module ElementTests using JuliaFEM.Test using JuliaFEM.Core: AbstractElement, Element, Field, FieldSet, test_element -using JuliaFEM.Core: Tri3 +using JuliaFEM.Core: Tri3, Quad4 import JuliaFEM.Core: get_basis, get_dbasis, calculate_normal_tangential_coordinates! import Base: size @@ -91,4 +91,15 @@ function test_calculate_normal_tangential_coordinates() end #test_calculate_normal_tangential_coordinates() +function test_manifold_determinant() + el = Quad4([1, 2, 3, 4]) + #el["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + el["geometry"] = Vector{Float64}[[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 1.0, 1.0], [0.0, 1.0, 1.0]] + # mother element area = 2*2 = 4, this element is 1, determinant should be 1/4 everywhere + d = det(el, [0.1, 0.2], 0.0) + d_expected = 0.25 + @test d == d_expected +end +#test_manifold_determinant() + end diff --git a/test/test_linear_elasticity.jl b/test/test_linear_elasticity.jl index b273fe7..aaa5e2c 100644 --- a/test/test_linear_elasticity.jl +++ b/test/test_linear_elasticity.jl @@ -36,13 +36,13 @@ function test_plane_stress_linear_elasticity_with_surface_load() free_dofs = Int64[3, 5, 6, 8] - info("initial force vector") ass = assemble(problem, 0.0) f = full(ass.force_vector) K = full(ass.stiffness_matrix) - dump(reshape(f, 2, 4)) - info("initial stiffness matrix") - dump(round(Int, K)[free_dofs, free_dofs]) +# info("initial force vector") +# dump(reshape(f, 2, 4)) +# info("initial stiffness matrix") +# dump(round(Int, K)[free_dofs, free_dofs]) u = zeros(2, 4) u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs] @@ -75,7 +75,7 @@ function test_continuum_elasticity_with_surface_load() end element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) set_geometry!(element1, nodes) - element1["youngs modulus"] = 9000.0 + element1["youngs modulus"] = 900.0 element1["poissons ratio"] = 0.25 element2 = Quad4([5, 6, 7, 8]) @@ -100,13 +100,13 @@ function test_continuum_elasticity_with_surface_load() free_dofs = find(vec(free_dofs')) info("free dofs: $free_dofs") - info("initial force vector") ass = assemble(problem, 0.0) f = full(ass.force_vector) K = full(ass.stiffness_matrix) - dump(reshape(f, 3, 8)) - info("initial stiffness matrix") - dump(round(Int, K)[free_dofs, free_dofs]) +# info("initial force vector") +# dump(reshape(f, 3, 8)) +# info("initial stiffness matrix") +# dump(round(Int, K)[free_dofs, free_dofs]) u = zeros(3, 8) u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs] diff --git a/test/test_potential_energy.jl b/test/test_potential_energy.jl index 0dbbe37..e311000 100644 --- a/test/test_potential_energy.jl +++ b/test/test_potential_energy.jl @@ -7,7 +7,7 @@ using JuliaFEM.Test using JuliaFEM.Core: AbstractProblem, Problem using JuliaFEM.Core: Element, Seg2, Quad4 -using JuliaFEM.Core: IntegrationPoint, solve! +using JuliaFEM.Core: IntegrationPoint, solve!, get_jacobian import JuliaFEM.Core: get_unknown_field_name, get_unknown_field_type, get_potential_energy @@ -34,7 +34,9 @@ function get_potential_energy(problem::Problem{HeatProblem}, element::Element{Qu gradT = element("temperature", ip, time, Val{:grad}, variation) Wint = (k + c*T) * 1/2*vecdot(gradT, gradT) Wext = f*T - return Wint - Wext + W = Wint - Wext + J = get_jacobian(element, ip, time) + return W*det(J) end function get_potential_energy(problem::Problem{HeatProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing) @@ -45,7 +47,8 @@ function get_potential_energy(problem::Problem{HeatProblem}, element::Element{Se Wint = 0.0 Wext = q0*T W = Wint - Wext - return W + J = get_jacobian(element, ip, time) + return W*norm(J) end function test_potential_energy_method() diff --git a/test/test_solver.jl b/test/test_solver.jl index 4c5ef8c..356207a 100644 --- a/test/test_solver.jl +++ b/test/test_solver.jl @@ -51,7 +51,7 @@ function test_linearsolver() info("Temperature at point X = $X is T = $T") @test isapprox(T, 100.0) end -#test_basic() +#test_linearsolver() function test_solvers() K = [ @@ -88,7 +88,8 @@ function test_solvers() @test isapprox(u1, expected) u2, la2 = solve(K, f, C, g, Val{:CHOLMOD}) @test isapprox(u2, expected) - include(Pkg.dir("JuliaFEM"*"/src/petsc.jl")) + # FIXME: how to dynamically include packages only if they are installed? + #include(Pkg.dir("JuliaFEM"*"/src/petsc.jl")) u3, la3 = solve(K, f, C, g, Val{:PETSc_GMRES}) @test isapprox(u3, expected) end