From 283621f4b1e1aedf703b89b8d76ca622b34ab26d Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 22 May 2016 00:22:17 +0300 Subject: [PATCH] rewrite elements --- src/JuliaFEM.jl | 65 +++++++++++++++++++++++++---- src/common.jl | 33 ++++++++++++++- src/core.jl | 71 -------------------------------- src/elasticity.jl | 44 +++++++++----------- src/elements.jl | 17 ++++---- src/fields.jl | 5 +++ src/integrate.jl | 75 +++++++++++++++++++++------------- src/lagrange_macro.jl | 48 ++++++++-------------- src/preprocess_aster_reader.jl | 2 - src/types.jl | 6 --- 10 files changed, 187 insertions(+), 179 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index d0e2017..f167ca9 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -6,25 +6,70 @@ This is JuliaFEM -- Finite Element Package """ module JuliaFEM -# include("common.jl") +importall Base +using ForwardDiff +using JLD +autodiffcache = ForwardDiffCache() +# export derivative, jacobian, hessian -""" JuliaFEM Core module. """ -module Core +include("common.jl") +typealias Node Vector{Float64} + +include("fields.jl") +export DCTI +#include("basis.jl") # interpolation of discrete fields +#include("symbolic.jl") # a thin symbolic layer for fields +#include("types.jl") # type definitions + +### ELEMENTS ### +include("elements.jl") # common element routines +export Element +include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro +export Quad4 + +#include("hierarchical.jl") # P-elements +#include("mortar_elements.jl") # Mortar elements +#include("equations.jl") + +include("integrate.jl") # default integration points for elements +include("sparse.jl") + +include("problems.jl") # common problem routines +export Problem + +include("elasticity.jl") # elasticity equations +export Elasticity + +include("dirichlet.jl") +include("heat.jl") +export assemble + +### ASSEMBLY + SOLVE ### +include("assembly.jl") +include("solver_utils.jl") +include("solvers.jl") + +### MORTAR STUFF ### +include("mortar.jl") # mortar projection + +include("abaqus_reader_old.jl") + +# rest of things +include("utils.jl") include("core.jl") -end module API include("api.jl") - # export .... - end module Preprocess +#= macro debug(msg) haskey(ENV, "DEBUG") || return return msg end +=# include("abaqus_reader.jl") include("preprocess_aster_reader.jl") end @@ -35,7 +80,13 @@ end """ JuliaFEM testing routines. """ module Test -include("test.jl") +if VERSION >= v"0.5-" + using Base.Test +else + using BaseTestNext +end +export @test, @testset, @test_throws +#include("test.jl") end module MaterialModels diff --git a/src/common.jl b/src/common.jl index 3f1c525..d78deb0 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,8 +1,6 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -# common routines - """ A very simple debugging macro. It executes commands if environment variable DEBUG is set. @@ -48,5 +46,36 @@ function set_debug_off!() pop!(ENV, "DEBUG"); end +""" Simple linspace extension to arrays. + +Examples +-------- +>>> linspace([0.0], [1.0], 3) +3-element Array{Array{Float64,1},1}: + [0.0] + [0.5] + [1.0] + +""" +function linspace{T<:Array}(X1::T, X2::T, n) + [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)] +end + +function resize!(A::SparseMatrixCSC, m::Int64, n::Int64) + (n == A.n) && (m == A.m) && return + @assert n >= A.n + @assert m >= A.m + append!(A.colptr, A.colptr[end]*ones(Int, m-A.m)) + A.n = n + A.m = m +end + +function ForwardDiff.derivative{T}(f::Function, S::Matrix{T}, args...) + shape = size(S) + wrapper(S::Vector) = f(reshape(S, shape)) + deriv = ForwardDiff.gradient(wrapper, vec(S), args...) + return reshape(deriv, shape) +end + export @debug, set_debug_on!, set_debug_off! diff --git a/src/core.jl b/src/core.jl index b90c9ab..eccdca1 100644 --- a/src/core.jl +++ b/src/core.jl @@ -1,73 +1,2 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -using JuliaFEM - -import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof, append! - -using ForwardDiff -autodiffcache = ForwardDiffCache() -# export derivative, jacobian, hessian - -using JLD - -""" Simple linspace extension to arrays. - -Examples --------- ->>> linspace([0.0], [1.0], 3) -3-element Array{Array{Float64,1},1}: - [0.0] - [0.5] - [1.0] - -""" -function Base.linspace{T<:Array}(X1::T, X2::T, n) - [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)] -end - -function Base.resize!(A::SparseMatrixCSC, m::Int64, n::Int64) - (n == A.n) && (m == A.m) && return - @assert n >= A.n - @assert m >= A.m - append!(A.colptr, A.colptr[end]*ones(Int, m-A.m)) - A.n = n - A.m = m -end - - -# fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb -include("fields.jl") -#include("basis.jl") # interpolation of discrete fields -include("symbolic.jl") # a thin symbolic layer for fields -include("types.jl") # type definitions - -### ELEMENTS ### -include("elements.jl") -include("lagrange.jl") # Lagrange elements -#include("hierarchical.jl") # P-elements -#include("mortar_elements.jl") # Mortar elements - -### EQUATIONS ### -include("integrate.jl") # default integration points for elements -include("sparse.jl") -include("problems.jl") -include("equations.jl") - -### FORMULATIION ### -include("dirichlet.jl") -include("heat.jl") -include("elasticity.jl") - -### ASSEMBLY + SOLVE ### -include("assembly.jl") -include("solver_utils.jl") -include("solvers.jl") - -### MORTAR STUFF ### -include("mortar.jl") # mortar projection - -include("abaqus_reader_old.jl") - -# rest of things -include("utils.jl") diff --git a/src/elasticity.jl b/src/elasticity.jl index b564aef..db58b48 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -43,29 +43,32 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El add!(assembly.f, gdofs, f) end +function assemble(problem::Problem{Elasticity}, element::Element, time=0.0) + 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::Real, ::Type{Val{:plane}}) +function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane_stress}}) props = problem.properties dim = get_unknown_field_dimension(problem) - nnodes = size(element, 2) + nnodes = length(element) BL = zeros(3, dim*nnodes) BNL = zeros(4, dim*nnodes) 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) - w = ip.weight*det(J) - N = element(ip, time) - dN = element(ip, time, Val{:grad}) + J = element(xi, time, Val{:Jacobian}) + w = w*det(J) + N = element(xi, time) + dN = element(xi, time, Val{:Grad}) # kinematics; calculate deformation gradient and strain gradu = zeros(dim, dim) if haskey(element, "displacement") - gradu += element("displacement", ip, time, Val{:grad}) + gradu += element("displacement", xi, time, Val{:Grad}) end strain = zeros(dim , dim) strain += 1/2*(gradu' + gradu) @@ -77,21 +80,12 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem # constitutive equations; material model (isotropic linear material here) # get_material(problem, element, ...) - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - 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 2d formulation: $(props.formulation)") - end + 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] # calculate stress S = D*[strain[1,1]; strain[2,2]; 2*strain[1,2]] @@ -127,8 +121,8 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem # volume load if haskey(element, "displacement load") - T = element("displacement load", ip, time) - f += vec(w*T*N) + b = element("displacement load", xi, time) + f += vec(w*N'*b) end end diff --git a/src/elements.jl b/src/elements.jl index bafa2d6..c25ee10 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -93,7 +93,16 @@ function update!(element::Element, field_name::ASCIIString, data::Union{Real, Ve element[field_name] = data 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)) + return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)' +end +""" Check existence of field. """ +function haskey(element::Element, field_name) + haskey(element.fields, field_name) +end #= @@ -205,10 +214,6 @@ function call(element::Element, field_name::ASCIIString, time::Number) return element[field_name](time) end -function get_dbasis{E<:AbstractElement}(::Type{E}, xi::Vector) - basis(xi) = vec(get_basis(E, xi)) - return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)' -end function get_basis{E}(element::Element{E}, ip::IntegrationPoint) return get_basis(E, ip.xi) @@ -339,10 +344,6 @@ function get_jacobian{E}(element::Element{E}, ip::IntegrationPoint, time::Real, end -""" Check does field exist. """ -function Base.haskey(element::Element, what) - haskey(element.fields, what) -end """ Calculate local normal-tangential coordinates for element. """ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real) diff --git a/src/fields.jl b/src/fields.jl index 6b8ee8a..3401b26 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -220,6 +220,11 @@ function Base.(:*){T<:Real}(c::T, field::DVTI) return DVTI(c*field.data) end +""" Interpolate with basis functions. """ +function Base.(:*)(N::Matrix, f::DCTI) + return f.data*N' +end + """ Multiply DVTI field with another vector T. Vector length must match to the field length and this can be used mainly for interpolation purposes, i.e., u = ∑ Nᵢuᵢ diff --git a/src/integrate.jl b/src/integrate.jl index e13612d..a891446 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -119,41 +119,60 @@ function get_integration_points(::Type{Tri3}) end -function get_integration_points(::Type{Quad4}, ::Type{Val{2}}) - [ - IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0), - IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0) - ] +function get_integration_points(::Type{Val{1}}) + return [2.0], [0.0] end -function get_integration_points(::Type{Quad4}, ::Type{Val{3}}) - p = [-sqrt(3/5), 0.0, sqrt(3/5)] - w = [5/9, 8/9, 5/9] - pts = vec([IntegrationPoint([p[i], p[j]], w[i]*w[j]) for i=1:3, j=1:3]) - return pts +function get_integration_points(::Type{Val{2}}) + return [1.0, 1.0], sqrt(1.0/3.0)*[-1.0, 1.0] end -function get_integration_points(::Type{Quad4}, ::Type{Val{5}}) - p = [ - -1/3*sqrt(5 + 2*sqrt(10/7)), - -1/3*sqrt(5 - 2*sqrt(10/7)), +function get_integration_points(::Type{Val{3}}) + return 1.0/9.0*[5.0, 8.0, 5.0], sqrt(3.0/5.0)*[-1.0, 0.0, 1.0] +end + +function get_integration_points(::Type{Val{4}}) + weights = 1.0/36.0*[ + 18.0+sqrt(30.0), + 18.0+sqrt(30.0), + 18.0-sqrt(30.0), + 18.0-sqrt(30.0)] + points = [ + sqrt( 3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)), + sqrt(-3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)), + sqrt( 3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)), + sqrt(-3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0))] + return weights, points +end + +function get_integration_points(::Type{Val{5}}) + weights = [ + 1.0/900.0*(322.0-13.0*sqrt(70.0)), + 1.0/900.0*(322.0+13.0*sqrt(70.0)), + 128.0/225.0, + 1.0/900.0*(322.0+13.0*sqrt(70.0)), + 1.0/900.0*(322.0-13.0*sqrt(70.0))] + points = [ + -1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0)), + -1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)), 0.0, - 1/3*sqrt(5 - 2*sqrt(10/7)), - 1/3*sqrt(5 + 2*sqrt(10/7))] - w = [ - (322-13*sqrt(70))/900, - (322+13*sqrt(70))/900, - 128/225, - (322+13*sqrt(70))/900, - (322-13*sqrt(70))/900] - pts = vec([IntegrationPoint([p[i], p[j]], w[i]*w[j]) for i=1:5, j=1:5]) - return pts + 1.0/3.0*sqrt(5.0 - 2.0*sqrt(10.0/7.0)), + 1.0/3.0*sqrt(5.0 + 2.0*sqrt(10.0/7.0))] + return weights, points end -function get_integration_points(::Type{Quad4}) - return get_integration_points(Quad4, Val{2}) +function get_integration_points(element::Quad4) + return get_integration_points(element, Val{2}) +end + +function get_integration_points(element::Quad4, ::Type{Val{2}}) + w, xi = get_integration_points(Val{2}) + [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:2, j=1:2 ] +end + +function get_integration_points(element::Quad4, ::Type{Val{3}}) + w, xi = get_integration_points(Val{3}) + [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:3, j=1:3 ] end ### 3d elements diff --git a/src/lagrange_macro.jl b/src/lagrange_macro.jl index af10ab0..cf75567 100644 --- a/src/lagrange_macro.jl +++ b/src/lagrange_macro.jl @@ -1,9 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -# Lagrange (Continous Galerkin) finite elements - -abstract CG <: AbstractElement +# Lagrange (Continous Galerkin) finite elements shape functions generated using macro. """ Given polynomial P and coordinates of reference element, calculate @@ -15,15 +13,6 @@ function calculate_lagrange_basis_coefficients(P, X) for i=1:nbasis A[i,:] = P(X[:, i]) end -# invA = inv(A)' -# basis(xi) = (invA*P(xi))' -# dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))' -# basis, dbasisdxi -# info(inv(A)) -# info(P([0.0, 0.0])) -# invA = inv(A)' -# basis(xi) = invA*P(xi) -# return basis return inv(A)' end @@ -41,21 +30,31 @@ Examples macro create_lagrange_element(element_name, element_description, X, P) eltype = esc(element_name) quote - global get_basis, get_dbasis, + global get_basis, length, size + #= get_reference_element_coordinates, get_reference_element_midpoint + =# + + type $eltype <: AbstractElement + end - #basis, dbasis = calculate_lagrange_basis($P, $X) C = calculate_lagrange_basis_coefficients($P, $X) basis(xi) = C*$P(xi) -# dbasis = ForwardDiff.jacobian(basis) - abstract $eltype <: CG - - function get_basis(::Type{$eltype}, xi::Vector) + function get_basis(element::$eltype, xi::Vector, time) return basis(xi)' end + function size(element::$eltype) + return size($X) + end + + function length(element::$eltype) + return size($X, 2) + end + + #= XX = refcoords($X) function get_reference_element_coordinates(::Type{$eltype}) return XX @@ -65,19 +64,12 @@ macro create_lagrange_element(element_name, element_description, X, P) function get_reference_element_midpoint(::Type{$eltype}) return XXX end -#= - function get_dbasis(::Type{$eltype}, xi::Vector{Float64}) - return dbasis(xi)' - end -=# function $eltype(args...) return Element{$eltype}(args...) end - function Base.size(::Type{$eltype}) - return Base.size($X) - end + =# end end @@ -97,10 +89,6 @@ end 0.0 0.0 1.0], (xi) -> [1.0, xi[1], xi[2]]) -#function get_reference_element_midpoint(::Type{Tri3}) -# return [1.0/3.0, 1.0/3.0] -#end - @create_lagrange_element(Tri6, "6 node quadratic triangle element", [0.0 1.0 0.0 0.5 0.5 0.0 0.0 0.0 1.0 0.0 0.5 0.5], diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 91ea7da..c1a2803 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -3,8 +3,6 @@ using HDF5 using JuliaFEM -using JuliaFEM.Core: Element, Quad4, Tri3, Tet4, Seg2, Hex8, update! - # TODO: this should be elsewhere function aster_create_elements(mesh, element_set, element_type=nothing; reverse_connectivity=false) diff --git a/src/types.jl b/src/types.jl index 05878ac..4b4683c 100644 --- a/src/types.jl +++ b/src/types.jl @@ -3,12 +3,6 @@ typealias Node Vector{Float64} -function ForwardDiff.derivative{T}(f::Function, S::Matrix{T}, args...) - shape = size(S) - wrapper(S::Vector) = f(reshape(S, shape)) - deriv = ForwardDiff.gradient(wrapper, vec(S), args...) - return reshape(deriv, shape) -end """ Integration point