From 53d640f477a8876604989bc1663589fe611b5fa8 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 22 May 2016 02:34:38 +0300 Subject: [PATCH] some tests pass now --- src/JuliaFEM.jl | 16 +- src/dirichlet.jl | 11 +- src/elasticity.jl | 47 ++-- src/elements.jl | 80 +++--- src/equations.jl | 13 - src/integrate.jl | 241 ++++++++---------- src/preprocess_aster_reader.jl | 3 +- src/problems.jl | 14 + ..._elasticity_2d_linear_with_surface_load.jl | 5 +- 9 files changed, 213 insertions(+), 217 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index f167ca9..8f0a144 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -23,9 +23,9 @@ export DCTI ### ELEMENTS ### include("elements.jl") # common element routines -export Element +export Element, update! include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro -export Quad4 +export Seg2, Tri3, Quad4, Hex8, Tet4 #include("hierarchical.jl") # P-elements #include("mortar_elements.jl") # Mortar elements @@ -41,13 +41,18 @@ include("elasticity.jl") # elasticity equations export Elasticity include("dirichlet.jl") +export Dirichlet + include("heat.jl") +export Heat + export assemble ### ASSEMBLY + SOLVE ### include("assembly.jl") include("solver_utils.jl") include("solvers.jl") +export Solver ### MORTAR STUFF ### include("mortar.jl") # mortar projection @@ -64,14 +69,9 @@ include("api.jl") end module Preprocess -#= -macro debug(msg) - haskey(ENV, "DEBUG") || return - return msg -end -=# include("abaqus_reader.jl") include("preprocess_aster_reader.jl") +export aster_create_elements, parse_aster_med_file end module Postprocess diff --git a/src/dirichlet.jl b/src/dirichlet.jl index cb7e218..6902244 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -43,26 +43,25 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele end # right hand side - for ip in get_integration_points(element, Val{3}) - w = ip.weight - J = get_jacobian(element, ip, time) + for (w, xi) in get_integration_points(element, Val{3}) + J = element(xi, time, Val{:Jacobian}) 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) + N = element(xi, time) for i=1:field_dim ldofs = gdofs[i:field_dim:end] if haskey(element, field_name*" $i") - g = element(field_name*" $i", ip, time) + g = element(field_name*" $i", xi, time) if get_formulation_type(problem) == :incremental # if having incremental formulation need to add previous # displacement to rhs (solving increment Δu ! haskey(element, "displacement") || continue - g_prev = element(field_name, ip, time) + g_prev = element(field_name, xi, time) g -= g_prev[i] end add!(assembly.g, ldofs, w*g*Ae*N') diff --git a/src/elasticity.jl b/src/elasticity.jl index db58b48..4557022 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -44,11 +44,15 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El end function assemble(problem::Problem{Elasticity}, element::Element, time=0.0) - assemble(problem, element, time, Val{problem.properties.formulation}) + problem.properties + if problem.properties.formulation in [:plane_stress, :plane_strain] + return assemble(problem, element, time, Val{:plane}) + end + return assemble(problem, element, time, Val{problem.properties.formulation}) end """ Elasticity equations for 2d cases. """ -function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane_stress}}) +function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) props = problem.properties dim = get_unknown_field_dimension(problem) @@ -82,10 +86,19 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem # get_material(problem, element, ...) E = element("youngs modulus", xi, time) nu = element("poissons ratio", xi, time) - D = E/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] + if props.formulation == :plane_stress + D = E/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + elseif props.formulation == :plane_strain + D = E/((1+nu)*(1-2*nu)) .* [ + 1-nu nu 0 + nu 1-nu 0 + 0 0 (1-2*nu)/2] + else + error("unknown plane formulation: $(props.formulation)") + end # calculate stress S = D*[strain[1,1]; strain[2,2]; 2*strain[1,2]] @@ -138,30 +151,30 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E Kt = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) - for ip in get_integration_points(element) + for (w, xi) in get_integration_points(element) - J = get_jacobian(element, ip, time) - N = element(ip, time) - w = ip.weight*norm(J) + J = element(xi, time, Val{:Jacobian}) + detJ = norm(J) + N = element(xi, time) if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - f += vec(w*T*N) + T = element("displacement traction force", xi, time) + f += w*vec(T*N)*detJ end for i=1:dim # traction force for ith component if haskey(element, "displacement traction force $i") - T = element("displacement traction force $i", ip, time) - f[i:dim:end] += vec(w*T*N) + T = element("displacement traction force $i", xi, time) + f[i:dim:end] += w*vec(T*N)*detJ end end if haskey(element, "nt displacement traction force") # traction force given in normal-tangential direction - T = element("nt displacement traction force", ip, time) - Q = element("normal-tangential coordinates", ip, time) - f += vec(w*Q'*T*N) + T = element("nt displacement traction force", xi, time) + Q = element("normal-tangential coordinates", xi, time) + f += w*vec(Q'*T*N)*detJ end end diff --git a/src/elements.jl b/src/elements.jl index c25ee10..8d486ae 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -74,6 +74,10 @@ function size{E}(element::Element{E}) size(element.properties) end +function size(element::Element, dim::Int) + size(element)[dim] +end + """ Update element field based on a dictionary of nodal data and connectivity information. Examples @@ -89,10 +93,16 @@ function update!(element::Element, field_name::ASCIIString, data::Dict) element[field_name] = [data[i] for i in get_connectivity(element)] end -function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair}...) +function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair}) element[field_name] = data end +function update!(elements::Vector, field_name::ASCIIString, data) + for element in elements + update!(element, field_name, data) + end +end + """ Evaluate partial derivatives of basis functions using ForwardDiff. """ function get_dbasis(element::Element, xi::Vector, time) basis(xi) = vec(get_basis(element, xi, time)) @@ -104,6 +114,36 @@ function haskey(element::Element, field_name) haskey(element.fields, field_name) end +function get_connectivity(element::Element) + return element.connectivity +end + +function get_gdofs(element::Element) + return get_gdofs(element, 1) +end + +""" Return dual basis transformation matrix Ae. """ +function get_dualbasis(element::Element, time) + nnodes = length(element) + De = zeros(nnodes, nnodes) + Me = zeros(nnodes, nnodes) + for (w, xi) in get_integration_points(element, Val{3}) + J = element(xi, time, Val{:Jacobian}) + JT = transpose(J) + if size(JT, 2) == 1 # plane problem + # || ∂X/∂ξ || + w *= norm(JT) + else + # || ∂X/∂ξ₁ × ∂X/∂ξ₂ || + w *= norm(cross(JT[:,1], JT[:,2])) + end + N = element(xi, time) + De += w*diagm(vec(N)) + Me += w*N'*N + end + return De, Me, De*inv(Me) +end + #= type Element{E} @@ -172,9 +212,6 @@ function Base.setindex!(element::Element, data::Tuple, name::ASCIIString) element.fields[name] = Field(data...) end -function get_connectivity(el::Element) - return el.connectivity -end typealias VecOrIP Union{Vector, IntegrationPoint} @@ -248,9 +285,6 @@ function find_elements(elements, nodes) return collect(s) end -function get_gdofs(element::Element) - return get_gdofs(element, 1) -end function get_dbasis{E}(element::Element{E}, ip::IntegrationPoint) return get_dbasis(E, ip.xi) @@ -260,33 +294,6 @@ function get_basis{E, T<:Real}(element::Element{E}, xi::T) return get_basis(E, xi) end -""" Return dual basis transformation matrix Ae. """ -function get_dualbasis(element::Element, time::Real) - if length(element.A) == 0 - nnodes = size(element, 2) - De = zeros(nnodes, nnodes) - Me = zeros(nnodes, nnodes) - for ip in get_integration_points(element, Val{3}) - w = ip.weight - J = get_jacobian(element, ip, time) - JT = transpose(J) - if size(JT, 2) == 1 # plane problem - # || ∂X/∂ξ || - w *= norm(JT) - else - # || ∂X/∂ξ₁ × ∂X/∂ξ₂ || - w *= norm(cross(JT[:,1], JT[:,2])) - end - N = element(ip, time) - De += w*diagm(vec(N)) - Me += w*N'*N - end - element.D = De - element.M = Me - element.A = De*inv(Me) - end - return element.D, element.M, element.A -end function call(element::Element, xi::VecOrIP, time::Real, ::Type{Val{:dualbasis}}) De, Me, Ae = get_dualbasis(element, time) @@ -486,10 +493,5 @@ function update!{T}(elements::Vector{Element{T}}, field_name::ASCIIString, data. update!(element, field_name, data...) end end -function update!(elements::Vector{Element}, field_name::ASCIIString, data...) - for element in elements - update!(element, field_name, data...) - end -end =# diff --git a/src/equations.jl b/src/equations.jl index 603e4e4..0851341 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -41,19 +41,6 @@ function has_residual_vector(problem::Problem, element::Element) return method_exists(get_residual_vector, default_args) end -function get_gdofs(element::Element, dim::Int) - conn = get_connectivity(element) - gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) - return gdofs -end - -function get_gdofs(element::Element, problem::Problem) - return get_gdofs(element, problem.dimension) -end - -function get_gdofs(problem::Problem, element::Element) - return get_gdofs(element, problem.dimension) -end """ Assemble element. """ function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number) diff --git a/src/integrate.jl b/src/integrate.jl index a891446..fd8579a 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -1,123 +1,9 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -# Let's drop here all integration schemes and some defaults for different element types -# maybe parse from txt file ..? - -### 1d elements - -typealias LineElement Union{Type{Seg2}, Type{Seg3}} - -function get_integration_points(::LineElement, ::Type{Val{1}}) - [ - IntegrationPoint([0.0], 2.0) - ] -end - -function get_integration_points(::LineElement, ::Type{Val{2}}) - [ - IntegrationPoint([-sqrt(1/3)], 1) - IntegrationPoint([+sqrt(1/3)], 1) - ] -end - -function get_integration_points(::LineElement, ::Type{Val{3}}) - [ - IntegrationPoint([0.0], 8/9), - IntegrationPoint([-sqrt(3/5)], 5/9), - IntegrationPoint([+sqrt(3/5)], 5/9) - ] -end - -function get_integration_points(::LineElement, ::Type{Val{4}}) - [ - IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36) - IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36) - IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36) - IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36) - ] -end - -function get_integration_points(::LineElement, ::Type{Val{5}}) - [ - IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900), - IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900), - IntegrationPoint([0.0], 128/225), - IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900), - IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900) - ] -end - -function get_integration_points(::Type{Seg2}) - return get_integration_points(Seg2, Val{2}) -end - -function get_integration_points(::Type{Seg3}) - return get_integration_points(Seg3, Val{3}) -end - -### 2d triangular elements - -# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF - -typealias TriangularElement Union{Type{Tri3}, Type{Tri6}} - -function get_integration_points(::TriangularElement, ::Type{Val{1}}) - # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html - [ - IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5) - ] -end - -function get_integration_points(::TriangularElement, ::Type{Val{2}}) - # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html - [ - IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0), - IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0), - IntegrationPoint([1.0/6.0, 1.0/6.0], 1.0/6.0) - ] -end - -function get_integration_points(::TriangularElement, ::Type{Val{3}}) - [ - IntegrationPoint([1/3, 1/3], 0.5*-0.5625), - IntegrationPoint([0.2, 0.2], 0.5*0.5208333333333333), - IntegrationPoint([0.2, 0.6], 0.5*0.5208333333333333), - IntegrationPoint([0.6, 0.2], 0.5*0.5208333333333333), - ] -end - -function get_integration_points(::TriangularElement, ::Type{Val{4}}) - # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF - # FIXME: something wrong here with weights ..? - [ - IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801), - IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801), - IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801), - IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532), - IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532), - IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532) - ] -end - -function get_integration_points(::TriangularElement, ::Type{Val{5}}) - # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF - # FIXME: something wrong here with weights ..? - [ - IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000), - IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851), - IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851), - IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851), - IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483), - IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483), - IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483) - ] -end - -function get_integration_points(::Type{Tri3}) - return get_integration_points(Tri3, Val{1}) -end +# Let's drop here all integration schemes and some defaults for different element types maybe parse from txt file ..? +### Gauss quadrature rules for one dimension function get_integration_points(::Type{Val{1}}) return [2.0], [0.0] @@ -161,8 +47,30 @@ function get_integration_points(::Type{Val{5}}) return weights, points end -function get_integration_points(element::Quad4) - return get_integration_points(element, Val{2}) +### "cartesian" elements, integration rules comes from tensor product + +### 1d elements + +typealias LineElement Union{Seg2, Seg3} + +function get_integration_points(element::LineElement, ::Type{Val{1}}) + w, xi = get_integration_points(Val{1}) + [ (w[i], [xi[i]]) for i=1:1 ] +end + +function get_integration_points(element::LineElement, ::Type{Val{2}}) + w, xi = get_integration_points(Val{2}) + [ (w[i], [xi[i]]) for i=1:2 ] +end + +function get_integration_points(element::LineElement, ::Type{Val{3}}) + w, xi = get_integration_points(Val{3}) + [ (w[i], [xi[i]]) for i=1:3 ] +end + +function get_integration_points{E<:LineElement}(element::Element{E}, ::Type{Val{3}}) + w, xi = get_integration_points(Val{3}) + [ (w[i], [xi[i]]) for i=1:3 ] end function get_integration_points(element::Quad4, ::Type{Val{2}}) @@ -175,19 +83,94 @@ function get_integration_points(element::Quad4, ::Type{Val{3}}) [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:3, j=1:3 ] end +function get_integration_points(element::Hex8, ::Type{Val{2}}) + w, xi = get_integration_points(Val{2}) + [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:2, j=1:2, k=1:2 ] +end + +### default number of integration points for each element + +function get_integration_points(element::Seg2) + get_integration_points(element, Val{2}) +end + +function get_integration_points(element::Seg3) + get_integration_points(element, Val{3}) +end + +function get_integration_points(element::Quad4) + get_integration_points(element, Val{2}) +end + +function get_integration_points(element::Hex8) + get_integration_points(element, Val{2}) +end + +### triangular and tetrahedral elements + +# http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + +typealias TriangularElement Union{Type{Tri3}, Type{Tri6}} + +function get_integration_points(::TriangularElement, ::Type{Val{1}}) + # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html + weights = [0.5] + points = Vector{Float64}[1.0/3.0*[1.0, 1.0]] + return weights, points +end + +function get_integration_points(::TriangularElement, ::Type{Val{2}}) + # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html + weights = 1.0/6.0*[1.0, 1.0, 1.0] + points = Vector{Float64}[ + [2.0/3.0, 1.0/6.0], + [1.0/6.0, 2.0/3.0], + [1.0/6.0, 1.0/6.0]] + return weights, points +end + +function get_integration_points(::TriangularElement, ::Type{Val{3}}) + [ + IntegrationPoint([1/3, 1/3], 0.5*-0.5625), + IntegrationPoint([0.2, 0.2], 0.5*0.5208333333333333), + IntegrationPoint([0.2, 0.6], 0.5*0.5208333333333333), + IntegrationPoint([0.6, 0.2], 0.5*0.5208333333333333), + ] +end + +function get_integration_points(::TriangularElement, ::Type{Val{4}}) + # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + # FIXME: something wrong here with weights ..? + [ + IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801), + IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801), + IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801), + IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532), + IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532), + IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532) + ] +end + +function get_integration_points(::TriangularElement, ::Type{Val{5}}) + # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + # FIXME: something wrong here with weights ..? + [ + IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000), + IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851), + IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851), + IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851), + IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483), + IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483), + IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483) + ] +end + +function get_integration_points(::Type{Tri3}) + return get_integration_points(Tri3, Val{1}) +end + ### 3d elements - -function get_integration_points(::Type{Hex8}, ::Type{Val{2}}) - p = 1.0/sqrt(3.0)*[-1.0, 1.0] - w = [1.0, 1.0] - return vec([IntegrationPoint([p[i], p[j], p[k]], w[i]*w[j]) for i=1:2, j=1:2, k=1:2]) -end - -function get_integration_points(::Type{Hex8}) - return get_integration_points(Hex8, Val{2}) -end - function get_integration_points(::Type{Tet4}) # http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html [ diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index c1a2803..047467b 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -4,7 +4,6 @@ using HDF5 using JuliaFEM -# TODO: this should be elsewhere function aster_create_elements(mesh, element_set, element_type=nothing; reverse_connectivity=false) elements = Element[] mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8, :TE4 => Tet4) @@ -25,7 +24,7 @@ function aster_create_elements(mesh, element_set, element_type=nothing; reverse_ if reverse_connectivity elcon = reverse(elcon) end - element = mapping[eltype](elcon) + element = Element(mapping[eltype], elcon) push!(elements, element) end update!(elements, "geometry", mesh["nodes"]) diff --git a/src/problems.jl b/src/problems.jl index 0c65b6a..b7688e8 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -267,6 +267,20 @@ function push!(problem::Problem, element) push!(problem.elements, element) end +function get_gdofs(element::Element, dim::Int) + conn = get_connectivity(element) + gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...)) + return gdofs +end + +function get_gdofs(element::Element, problem::Problem) + return get_gdofs(element, problem.dimension) +end + +function get_gdofs(problem::Problem, element::Element) + return get_gdofs(element, problem.dimension) +end + """ Find dofs corresponding to nodes. """ function find_dofs_by_nodes(problem::Problem, nodes) dim = get_unknown_field_dimension(problem) diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl index 9779baf..ce6bb09 100644 --- a/test/test_elasticity_2d_linear_with_surface_load.jl +++ b/test/test_elasticity_2d_linear_with_surface_load.jl @@ -1,11 +1,10 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +using JuliaFEM +using JuliaFEM.Preprocess using JuliaFEM.Test -using JuliaFEM.Preprocess: aster_create_elements, parse_aster_med_file -using JuliaFEM.Core: Problem, Elasticity, Dirichlet, Solver, update! - @testset "test 2d linear elasticity with surface load" begin meshfile = "/geometry/2d_block/BLOCK_1elem.med" mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile)