diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 0b3cd47..e8906a0 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -6,13 +6,80 @@ abstract StandardBasis abstract DualBasis global const BiorthogonalBasis = DualBasis -function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]; basis=StandardBasis) - return BoundaryProblem{DirichletProblem{basis}}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements) +type Dirichlet <: AbstractProblem + dual_basis :: Bool end -function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]; basis=StandardBasis) - return BoundaryProblem{DirichletProblem{basis}}(problem_name, parent_field_name, parent_field_dim, dim, elements) + +function Dirichlet() + Dirichlet(true) end +function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichlet}, element::Element, time::Real) + + @assert problem.properties.dual_basis + + # get dimension and name of PARENT field + field_dim = problem.parent_field_dim + field_name = problem.parent_field_name + gdofs = get_gdofs(element, field_dim) + + # calculate bi-orthogonal basis transformation matrix Ae + nnodes = size(element, 2) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for ip in get_integration_points(element, Val{2}) + w = ip.weight + J = get_jacobian(element, ip, time) + JT = transpose(J) + if size(JT, 2) == 1 # plane problem + w *= norm(JT) + else + w *= norm(cross(JT[:,1], JT[:,2])) + end + N = element(ip, time) + De += w*diagm(vec(N)) + Me += w*N'*N + end + Ae = De*inv(Me) + + # do the actual integration + for ip in get_integration_points(element, Val{2}) + w = ip.weight + J = get_jacobian(element, ip, time) + JT = transpose(J) + if size(JT, 2) == 1 # plane problem + w *= norm(JT) + else + w *= norm(cross(JT[:,1], JT[:,2])) + end + N = element(ip, time) + Phi = (Ae*N')' + A = w*Phi'*N + A[abs(A) .< 1.0e-12] = 0 + + if haskey(element, field_name) + for i=1:field_dim + g = element(field_name, ip, time) + ldofs = gdofs[i:field_dim:end] + add!(assembly.C1, ldofs, ldofs, A) + add!(assembly.C2, ldofs, ldofs, A) + add!(assembly.g, ldofs, w*g*Phi') + end + else + for i=1:field_dim + ldofs = gdofs[i:field_dim:end] + if haskey(element, field_name*" $i") + g = element(field_name*" $i", ip, time) + add!(assembly.C1, ldofs, ldofs, A) + add!(assembly.C2, ldofs, ldofs, A) + add!(assembly.g, ldofs, w*g*Phi') + end + end + end + end +end + + function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Real) # get dimension and name of PARENT field @@ -61,74 +128,3 @@ function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichle end end -function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem{BiorthogonalBasis}}, element::Element, time::Real) - - # get dimension and name of PARENT field - field_dim = problem.parent_field_dim - field_name = problem.parent_field_name - gdofs = get_gdofs(element, field_dim) - - # calculate bi-orthogonal basis transformation matrix Ae - nnodes = size(element, 2) - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(element, Val{2}) - w = ip.weight - J = get_jacobian(element, ip, time) - JT = transpose(J) - if size(JT, 2) == 1 # plane problem - w *= norm(JT) - else - w *= norm(cross(JT[:,1], JT[:,2])) - end - N = element(ip, time) - De += w*diagm(vec(N)) - Me += w*N'*N - end - Ae = De*inv(Me) - - # do the actual integration - for ip in get_integration_points(element, Val{2}) - w = ip.weight - J = get_jacobian(element, ip, time) - JT = transpose(J) - if size(JT, 2) == 1 # plane problem - w *= norm(JT) - else - w *= norm(cross(JT[:,1], JT[:,2])) - end - N = element(ip, time) - Phi = (Ae*N')' - A = w*Phi'*N - A[abs(A) .< 1.0e-9] = 0 - - # C1 matrix is always the same -# for i=1:field_dim -# ldofs = gdofs[i:field_dim:end] -# add!(assembly.C1, ldofs, ldofs, A) -# end - - if haskey(element, field_name) - # add all dimensions at once if defined element["blaa"] = 0.0 - for i=1:field_dim - g = element(field_name, ip, time) - ldofs = gdofs[i:field_dim:end] - add!(assembly.C1, ldofs, ldofs, A) - add!(assembly.C2, ldofs, ldofs, A) - add!(assembly.g, ldofs, w*g*Phi') - end - else - for i=1:field_dim - ldofs = gdofs[i:field_dim:end] - if haskey(element, field_name*" $i") - g = element(field_name*" $i", ip, time) - add!(assembly.C1, ldofs, ldofs, A) - add!(assembly.C2, ldofs, ldofs, A) - add!(assembly.g, ldofs, w*g*Phi') -# else -# add!(assembly.D, ldofs, ldofs, A) - end - end - end - end -end diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl index 1201db6..c68eadd 100644 --- a/src/linear_elasticity.jl +++ b/src/linear_elasticity.jl @@ -3,6 +3,88 @@ # Linear elasticity +""" Concrete Elasticity type. """ +type Elasticity <: AbstractProblem + plane_stress :: Bool + nonlinear_geometry :: Bool +end +function Elasticity() + Elasticity(false, false) +end + +function get_unknown_field_name(::Type{Elasticity}) + return "displacement" +end + +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real) + + # assemble plane stress problem + if problem.properties.plane_stress + return assemble!(assembly, problem, element, time, Val{:plane_stress}) + end + +end + +""" Elasticity equations, plane stress. """ +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:plane_stress}}) + + gdofs = get_gdofs(element, problem.dim) + ndim, nnodes = size(element) + B = zeros(3, 2*nnodes) + for ip in get_integration_points(element) + 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) + E_ = element("youngs modulus", ip, time) + C = E_/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + dN = element(ip, time, Val{:grad}) + fill!(B, 0.0) + for i=1:size(dN, 2) + B[1, 2*(i-1)+1] = dN[1,i] + B[2, 2*(i-1)+2] = dN[2,i] + B[3, 2*(i-1)+1] = dN[2,i] + B[3, 2*(i-1)+2] = dN[1,i] + end + Kt = w*B'*C*B*det(J) + add!(assembly.stiffness_matrix, gdofs, gdofs, Kt) + end + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + 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*norm(J) + add!(assembly.force_vector, gdofs, vec(L)) + end + for dim in 1:problem.dim + if haskey(element, "displacement traction force $dim") + T = element("displacement traction force $dim", ip, time) + ldofs = gdofs[dim:problem.dim:end] + L = w*T*N*norm(J) + add!(assembly.force_vector, ldofs, vec(L)) + end + end + if haskey(element, "displacement traction force N") + # surface pressure + p = zeros(2) + p[1] = element("displacement traction force N", ip, time) + R = element("normal-tangential coordinates", ip, time) + T = R'*p + L = w*T*N*norm(J) + add!(assembly.force_vector, gdofs, vec(L)) + end + end +end + + + + abstract LinearElasticityProblem <: ElasticityProblem function LinearElasticityProblem(name="linear elasticity", dim::Int=3, elements=[]) @@ -84,76 +166,6 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble end end -abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem - -function PlaneStressLinearElasticityProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) - return Problem{PlaneStressLinearElasticityProblem}(name, dim, elements) -end - -""" Elasticity equations, plane stress. """ -function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) - - gdofs = get_gdofs(element, problem.dim) - ndim, nnodes = size(E) - B = zeros(3, 2*nnodes) - for ip in get_integration_points(element) - 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) - E_ = element("youngs modulus", ip, time) - C = E_/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 2*(i-1)+1] = dN[1,i] - B[2, 2*(i-1)+2] = dN[2,i] - B[3, 2*(i-1)+1] = dN[2,i] - B[3, 2*(i-1)+2] = dN[1,i] - end - Kt = w*B'*C*B*det(J) - add!(assembly.stiffness_matrix, gdofs, gdofs, Kt) - # solve residual, i.e. K du = K \ -(Ku(prev) - F) - # in first iteration u(prev) typically 0 -> no effect - # but if geometrical or material nonlinearities iterations are needed -# if haskey(element, "displacement") -# u_prev = element("displacement", time) -# u_prev = vec(u_prev) -# add!(assembly.force_vector, gdofs, -Kt*u_prev) -# end - end - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - 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*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - for dim in 1:problem.dim - if haskey(element, "displacement traction force $dim") - T = element("displacement traction force $dim", ip, time) - ldofs = gdofs[dim:problem.dim:end] - L = w*T*N*norm(J) - add!(assembly.force_vector, ldofs, vec(L)) - end - end - if haskey(element, "displacement traction force N") - # surface pressure - p = zeros(2) - p[1] = element("displacement traction force N", ip, time) - R = element("normal-tangential coordinates", ip, time) - T = R'*p - L = w*T*N*norm(J) - add!(assembly.force_vector, gdofs, vec(L)) - end - end -end ############################### # Plastic material # diff --git a/src/problems.jl b/src/problems.jl index ead4f30..3d228c2 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -50,10 +50,11 @@ type FieldProblem{T} dim :: Int elements :: Vector{Element} assembly :: FieldAssembly + properties :: T end -function FieldProblem(problem_type::DataType, name::ASCIIString, dim::Int, +function FieldProblem(problem::DataType, name::ASCIIString, dim::Int, elements=[]) - FieldProblem{problem_type}(name, dim, elements, FieldAssembly()) + FieldProblem{problem}(name, dim, elements, FieldAssembly(), problem()) end @@ -124,14 +125,15 @@ type BoundaryProblem{T} parent_field_dim :: Int elements :: Vector{Element} assembly :: BoundaryAssembly + properties :: T end -function BoundaryProblem(problem_type::DataType, +function BoundaryProblem(problem::DataType, name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, elements=[]) - BoundaryProblem{problem_type}(name, parent_field_name, parent_field_dim, - elements, BoundaryAssembly()) + BoundaryProblem{problem}(name, parent_field_name, parent_field_dim, + elements, BoundaryAssembly(), problem()) end function update!{P}(problem::BoundaryProblem{P}, solution::Vector{Float64}) diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl index 6b9530c..907f495 100644 --- a/test/test_elasticity_surface_load.jl +++ b/test/test_elasticity_surface_load.jl @@ -2,9 +2,8 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM.Test -using JuliaFEM.Core: Node, update!, Quad4, Seg2, PlaneStressLinearElasticityProblem, - assemble, FieldProblem, BoundaryProblem, DirichletProblem, - DirectSolver, DualBasis +using JuliaFEM.Core: Node, update!, Quad4, Seg2, assemble, BoundaryProblem, + Problem, Elasticity, Solver, Dirichlet @testset "test 2d linear elasticity with surface load." begin @@ -27,7 +26,8 @@ using JuliaFEM.Core: Node, update!, Quad4, Seg2, PlaneStressLinearElasticityProb update!(element1, "poissons ratio", nu) # update!(element2, "displacement traction force", [0.0, f]) update!(element2, "displacement traction force", Vector{Float64}[[0.0, f], [0.0, f]]) - elasticity_problem = FieldProblem(PlaneStressLinearElasticityProblem, "block", 2) + elasticity_problem = Problem(Elasticity, "block", 2) + elasticity_problem.properties.plane_stress = true push!(elasticity_problem, element1, element2) # dirichlet boundary condition, symmetry @@ -36,37 +36,15 @@ using JuliaFEM.Core: Node, update!, Quad4, Seg2, PlaneStressLinearElasticityProb update!([sym13, sym23], "geometry", nodes) update!(sym13, "displacement 2", 0.0) update!(sym23, "displacement 1", 0.0) - #sym23["displacement 1"] = 0.0 - #sym13["displacement 2"] = 0.0 # name, unknown field, unknown field dimension - boundary_problem = BoundaryProblem(DirichletProblem, "symmetry boundaries", "displacement", 2) + boundary_problem = BoundaryProblem(Dirichlet, "symmetry boundaries", "displacement", 2) push!(boundary_problem, sym13, sym23) - solver = DirectSolver("solve block problem") - solver.solve_residual = false + solver = Solver("solve block problem") push!(solver, elasticity_problem) push!(solver, boundary_problem) - #= - add_linear_system_solver_posthook!(solver, - function (solver) - info("solution vector") - dump(solver.x) - end) - =# - - call(solver, 0.0) - - #= - free_dofs = Int64[3, 5, 6, 8] - ass = assemble(problem, 0.0) - f = full(ass.force_vector) - K = full(ass.stiffness_matrix) - u = zeros(2, 4) - u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs] - info("result vector") - dump(u) - =# + call(solver) u_disp = element1("displacement", [1.0, 1.0], 0.0) info("Displacement = $u_disp")