mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Set up entry point for element assembly
From here we can preallocate matrices to improve efficiency of code.
This commit is contained in:
Jukka Aho
+243
-194
@@ -53,17 +53,19 @@ function get_formulation_type(problem::Problem{Elasticity})
|
||||
return :incremental
|
||||
end
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time=0.0)
|
||||
props = problem.properties
|
||||
gdofs = get_gdofs(problem, element)
|
||||
formulation = props.formulation
|
||||
if formulation in [:plane_stress, :plane_strain]
|
||||
formulation = :plane
|
||||
end
|
||||
Km, Kg, f = assemble(problem, element, time, Val{formulation})
|
||||
add!(assembly.K, gdofs, gdofs, Km)
|
||||
add!(assembly.Kg, gdofs, gdofs, Kg)
|
||||
add!(assembly.f, gdofs, f)
|
||||
"""
|
||||
Start finite element assembly procedure for Elasticity problem.
|
||||
"""
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time)
|
||||
assemble!(assembly, problem, elements, time, Val{problem.properties.formulation})
|
||||
end
|
||||
|
||||
"""
|
||||
This is for backward compatibility, will be removed asap.
|
||||
"""
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time)
|
||||
warn("try to avoid single element assembly function as it's not possible to preallocate causing a slow code")
|
||||
assemble!(assembly, problem, [element], time, Val{problem.properties.formulation})
|
||||
end
|
||||
|
||||
include("problems_elasticity_2d.jl")
|
||||
@@ -100,235 +102,282 @@ function get_keys(element)
|
||||
map(x -> all_keys[x], idx)
|
||||
end
|
||||
|
||||
""" Elasticity equations, 3d nonlinear. """
|
||||
function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
|
||||
""" Continuum elements assembly entry point.
|
||||
|
||||
This splits elements to arrays by their type and assemble one element type
|
||||
at time. This makes it possible to pre-allocate matrices common to same type
|
||||
of elements.
|
||||
"""
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity},
|
||||
all_elements::Vector{Element}, time, ::Type{Val{:continuum}})
|
||||
element_types = unique(map(get_element_type, all_elements))
|
||||
for element_type in element_types
|
||||
elements = filter_by_element_type(element_type, all_elements)
|
||||
# FIXME: there must be better way to do this
|
||||
# to promote array for certain elemene type
|
||||
elements = [element for element in elements]
|
||||
nelements = length(elements)
|
||||
debug("elasticity 3d: assembling $nelements of type $element_type")
|
||||
assemble!(assembly, problem, elements, time, Val{:continuum})
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
""" Assemble 3d continuum elements in general solid mechanics problem. """
|
||||
function assemble!{El<:Elasticity3DVolumeElements}(assembly::Assembly,
|
||||
problem::Problem{Elasticity},
|
||||
elements::Vector{Element{El}},
|
||||
time, ::Type{Val{:continuum}})
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = length(element)
|
||||
ndofs = dim*nnodes
|
||||
BL = zeros(6, ndofs)
|
||||
BNL = zeros(9, ndofs)
|
||||
Km = zeros(ndofs, ndofs)
|
||||
Kg = zeros(ndofs, ndofs)
|
||||
f = zeros(ndofs)
|
||||
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
for element in elements
|
||||
|
||||
# kinematics; calculate deformation gradient and strain
|
||||
nnodes = length(element)
|
||||
ndofs = dim*nnodes
|
||||
BL = zeros(6, ndofs)
|
||||
BNL = zeros(9, ndofs)
|
||||
Km = zeros(ndofs, ndofs)
|
||||
Kg = zeros(ndofs, ndofs)
|
||||
f = zeros(ndofs)
|
||||
|
||||
gradu = zeros(dim, dim)
|
||||
if haskey(element, "displacement")
|
||||
gradu += element("displacement", ip, time, Val{:Grad})
|
||||
end
|
||||
strain = 1/2*(gradu' + gradu)
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
|
||||
F = eye(dim)
|
||||
if props.finite_strain
|
||||
F += gradu
|
||||
strain += 1/2*gradu'*gradu
|
||||
end
|
||||
# kinematics; calculate deformation gradient and strain
|
||||
|
||||
# material stiffness start
|
||||
gradu = zeros(dim, dim)
|
||||
if haskey(element, "displacement")
|
||||
gradu += element("displacement", ip, time, Val{:Grad})
|
||||
end
|
||||
strain = 1/2*(gradu' + gradu)
|
||||
|
||||
fill!(BL, 0.0)
|
||||
for i=1:nnodes
|
||||
BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
|
||||
BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
|
||||
BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
|
||||
BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i]
|
||||
BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i]
|
||||
BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i]
|
||||
BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i]
|
||||
BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i]
|
||||
BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i]
|
||||
BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
|
||||
BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
|
||||
BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i]
|
||||
BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i]
|
||||
BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i]
|
||||
BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i]
|
||||
BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i]
|
||||
BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
|
||||
BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
|
||||
end
|
||||
F = eye(dim)
|
||||
if props.finite_strain
|
||||
F += gradu
|
||||
strain += 1/2*gradu'*gradu
|
||||
end
|
||||
|
||||
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
|
||||
# material stiffness start
|
||||
|
||||
# calculate stress
|
||||
E = element("youngs modulus", ip, time)
|
||||
nu = element("poissons ratio", ip, time)
|
||||
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
|
||||
1.0-nu nu nu 0.0 0.0 0.0
|
||||
nu 1.0-nu nu 0.0 0.0 0.0
|
||||
nu nu 1.0-nu 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.5-nu 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.5-nu 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.5-nu]
|
||||
fill!(BL, 0.0)
|
||||
for i=1:nnodes
|
||||
BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
|
||||
BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
|
||||
BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
|
||||
BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i]
|
||||
BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i]
|
||||
BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i]
|
||||
BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i]
|
||||
BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i]
|
||||
BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i]
|
||||
BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
|
||||
BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
|
||||
BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i]
|
||||
BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i]
|
||||
BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i]
|
||||
BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i]
|
||||
BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i]
|
||||
BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
|
||||
BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
|
||||
end
|
||||
|
||||
element_keys = get_keys(element)
|
||||
strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
|
||||
|
||||
if "plasticity" in element_keys
|
||||
plastic_def = element("plasticity")[ip.id]
|
||||
# calculate stress
|
||||
E = element("youngs modulus", ip, time)
|
||||
nu = element("poissons ratio", ip, time)
|
||||
D = E/((1.0+nu)*(1.0-2.0*nu)) * [
|
||||
1.0-nu nu nu 0.0 0.0 0.0
|
||||
nu 1.0-nu nu 0.0 0.0 0.0
|
||||
nu nu 1.0-nu 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.5-nu 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.5-nu 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.5-nu]
|
||||
|
||||
calculate_stress! = plastic_def["type"]
|
||||
yield_surface_ = plastic_def["yield_surface"]
|
||||
params = plastic_def["params"]
|
||||
element_keys = get_keys(element)
|
||||
|
||||
initialize_internal_params!(params, ip, Val{:type_3d})
|
||||
if "plasticity" in element_keys
|
||||
plastic_def = element("plasticity")[ip.id]
|
||||
|
||||
if time == 0.0
|
||||
error("Given step time = $(time). Please select time > 0.0")
|
||||
calculate_stress! = plastic_def["type"]
|
||||
yield_surface_ = plastic_def["yield_surface"]
|
||||
params = plastic_def["params"]
|
||||
|
||||
initialize_internal_params!(params, ip, Val{:type_3d})
|
||||
|
||||
if time == 0.0
|
||||
error("Given step time = $(time). Please select time > 0.0")
|
||||
end
|
||||
|
||||
t_last = ip("prev_time", time)
|
||||
update!(ip, "prev_time", time => t_last)
|
||||
|
||||
dt = time - t_last
|
||||
|
||||
stress_last = ip("stress", t_last)
|
||||
strain_last = ip("strain", t_last)
|
||||
|
||||
dstrain_vec = strain_vec - strain_last
|
||||
stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
Dtan = [0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0]
|
||||
calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})
|
||||
|
||||
|
||||
# plastic_def = element.dev["plasticity"]
|
||||
# calculate_stress! = plastic_def["stress"]
|
||||
# params = plastic_def["params"]
|
||||
# yield_surface_ = plastic_def["yield_surface"]
|
||||
# (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d})
|
||||
# dstrain_vec = strain_vec - strain_last
|
||||
|
||||
# calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d})
|
||||
else
|
||||
|
||||
stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
|
||||
Dtan = D
|
||||
end
|
||||
|
||||
:strain in props.store_fields && update!(ip, "strain", time => strain_vec)
|
||||
:stress in props.store_fields && update!(ip, "stress", time => stress_vec)
|
||||
:stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
|
||||
:stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
|
||||
:stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3])
|
||||
:stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
|
||||
:stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
|
||||
:stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
|
||||
:plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
|
||||
|
||||
Km += w*BL'*Dtan*BL
|
||||
# material stiffness end
|
||||
|
||||
if props.geometric_stiffness
|
||||
# take geometric stiffness into account
|
||||
|
||||
fill!(BNL, 0.0)
|
||||
|
||||
for i=1:size(dN, 2)
|
||||
BNL[1, 3*(i-1)+1] = dN[1,i]
|
||||
BNL[2, 3*(i-1)+1] = dN[2,i]
|
||||
BNL[3, 3*(i-1)+1] = dN[3,i]
|
||||
BNL[4, 3*(i-1)+2] = dN[1,i]
|
||||
BNL[5, 3*(i-1)+2] = dN[2,i]
|
||||
BNL[6, 3*(i-1)+2] = dN[3,i]
|
||||
BNL[7, 3*(i-1)+3] = dN[1,i]
|
||||
BNL[8, 3*(i-1)+3] = dN[2,i]
|
||||
BNL[9, 3*(i-1)+3] = dN[3,i]
|
||||
end
|
||||
|
||||
t_last = ip("prev_time", time)
|
||||
update!(ip, "prev_time", time => t_last)
|
||||
S3 = zeros(3*dim, 3*dim)
|
||||
S3[1,1] = stress_vec[1]
|
||||
S3[2,2] = stress_vec[2]
|
||||
S3[3,3] = stress_vec[3]
|
||||
S3[1,2] = S3[2,1] = stress_vec[4]
|
||||
S3[2,3] = S3[3,2] = stress_vec[5]
|
||||
S3[1,3] = S3[3,1] = stress_vec[6]
|
||||
S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
|
||||
|
||||
dt = time - t_last
|
||||
Kg += w*BNL'*S3*BNL
|
||||
|
||||
stress_last = ip("stress", t_last)
|
||||
strain_last = ip("strain", t_last)
|
||||
end
|
||||
|
||||
dstrain_vec = strain_vec - strain_last
|
||||
stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
Dtan = [0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0;
|
||||
0.0 0.0 0.0 0.0 0.0 0.0]
|
||||
calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})
|
||||
# external load start
|
||||
|
||||
if haskey(element, "displacement load")
|
||||
T = element("displacement load", ip, time)
|
||||
f += w*vec(T*N)
|
||||
end
|
||||
|
||||
# plastic_def = element.dev["plasticity"]
|
||||
# calculate_stress! = plastic_def["stress"]
|
||||
# params = plastic_def["params"]
|
||||
# yield_surface_ = plastic_def["yield_surface"]
|
||||
# (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d})
|
||||
# dstrain_vec = strain_vec - strain_last
|
||||
for i=1:dim
|
||||
if haskey(element, "displacement load $i")
|
||||
b = element("displacement load $i", ip, time)
|
||||
f[i:dim:end] += w*vec(b*N)
|
||||
end
|
||||
end
|
||||
|
||||
# calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d})
|
||||
else
|
||||
# external load end
|
||||
|
||||
if get_formulation_type(problem) == :incremental
|
||||
f -= w*BL'*stress_vec
|
||||
end
|
||||
|
||||
stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
|
||||
Dtan = D
|
||||
end
|
||||
|
||||
:strain in props.store_fields && update!(ip, "strain", time => strain_vec)
|
||||
:stress in props.store_fields && update!(ip, "stress", time => stress_vec)
|
||||
:stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
|
||||
:stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
|
||||
:stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3])
|
||||
:stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
|
||||
:stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
|
||||
:stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
|
||||
:plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
|
||||
gdofs = get_gdofs(problem, element)
|
||||
|
||||
Km += w*BL'*Dtan*BL
|
||||
# material stiffness end
|
||||
# add contributions to K, Kg, f
|
||||
add!(assembly.K, gdofs, gdofs, Km)
|
||||
|
||||
if props.geometric_stiffness
|
||||
# take geometric stiffness into account
|
||||
|
||||
fill!(BNL, 0.0)
|
||||
|
||||
for i=1:size(dN, 2)
|
||||
BNL[1, 3*(i-1)+1] = dN[1,i]
|
||||
BNL[2, 3*(i-1)+1] = dN[2,i]
|
||||
BNL[3, 3*(i-1)+1] = dN[3,i]
|
||||
BNL[4, 3*(i-1)+2] = dN[1,i]
|
||||
BNL[5, 3*(i-1)+2] = dN[2,i]
|
||||
BNL[6, 3*(i-1)+2] = dN[3,i]
|
||||
BNL[7, 3*(i-1)+3] = dN[1,i]
|
||||
BNL[8, 3*(i-1)+3] = dN[2,i]
|
||||
BNL[9, 3*(i-1)+3] = dN[3,i]
|
||||
end
|
||||
|
||||
S3 = zeros(3*dim, 3*dim)
|
||||
S3[1,1] = stress_vec[1]
|
||||
S3[2,2] = stress_vec[2]
|
||||
S3[3,3] = stress_vec[3]
|
||||
S3[1,2] = S3[2,1] = stress_vec[4]
|
||||
S3[2,3] = S3[3,2] = stress_vec[5]
|
||||
S3[1,3] = S3[3,1] = stress_vec[6]
|
||||
S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
|
||||
|
||||
Kg += w*BNL'*S3*BNL
|
||||
|
||||
add!(assembly.Kg, gdofs, gdofs, Kg)
|
||||
end
|
||||
|
||||
# external load start
|
||||
|
||||
if haskey(element, "displacement load")
|
||||
T = element("displacement load", ip, time)
|
||||
f += w*vec(T*N)
|
||||
end
|
||||
|
||||
for i=1:dim
|
||||
if haskey(element, "displacement load $i")
|
||||
b = element("displacement load $i", ip, time)
|
||||
f[i:dim:end] += w*vec(b*N)
|
||||
end
|
||||
end
|
||||
|
||||
# external load end
|
||||
|
||||
if get_formulation_type(problem) == :incremental
|
||||
f -= w*BL'*stress_vec
|
||||
end
|
||||
add!(assembly.f, gdofs, f)
|
||||
|
||||
end
|
||||
|
||||
return Km, Kg, f
|
||||
return nothing
|
||||
end
|
||||
|
||||
""" Elasticity equations, surface traction for continuum formulation. """
|
||||
function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
|
||||
function assemble!{El<:Elasticity3DSurfaceElements}(assembly::Assembly,
|
||||
problem::Problem{Elasticity},
|
||||
elements::Vector{Element{El}},
|
||||
time, ::Type{Val{:continuum}})
|
||||
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
nnodes = size(element, 2)
|
||||
Km = zeros(dim*nnodes, dim*nnodes)
|
||||
Kg = zeros(dim*nnodes, dim*nnodes)
|
||||
f = zeros(dim*nnodes)
|
||||
|
||||
has_concentrated_forces = false
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
if haskey(element, "displacement traction force")
|
||||
T = element("displacement traction force", ip, time)
|
||||
f += w*vec(T*N)
|
||||
end
|
||||
for i in 1:dim
|
||||
if haskey(element, "displacement traction force $i")
|
||||
T = element("displacement traction force $i", ip, time)
|
||||
f[i:dim:end] += w*vec(T*N)
|
||||
for element in elements
|
||||
nnodes = size(element, 2)
|
||||
f = zeros(dim*nnodes)
|
||||
|
||||
has_concentrated_forces = false
|
||||
for ip in get_integration_points(element)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
if haskey(element, "displacement traction force")
|
||||
T = element("displacement traction force", ip, time)
|
||||
f += w*vec(T*N)
|
||||
end
|
||||
if haskey(element, "concentrated force $i")
|
||||
has_concentrated_forces = true
|
||||
T = element("concentrated force $i", ip, time)
|
||||
f[i:dim:end] += w*vec(T*N)
|
||||
for i in 1:dim
|
||||
if haskey(element, "displacement traction force $i")
|
||||
T = element("displacement traction force $i", ip, time)
|
||||
f[i:dim:end] += w*vec(T*N)
|
||||
end
|
||||
if haskey(element, "concentrated force $i")
|
||||
has_concentrated_forces = true
|
||||
T = element("concentrated force $i", ip, time)
|
||||
f[i:dim:end] += w*vec(T*N)
|
||||
end
|
||||
end
|
||||
if haskey(element, "surface pressure")
|
||||
J = element(ip, time, Val{:Jacobian})'
|
||||
n = cross(J[:,1], J[:,2])
|
||||
n /= norm(n)
|
||||
# sign convention, positive pressure is towards surface
|
||||
p = -element("surface pressure", ip, time)
|
||||
f += w*p*vec(n*N)
|
||||
end
|
||||
end
|
||||
if haskey(element, "surface pressure")
|
||||
J = element(ip, time, Val{:Jacobian})'
|
||||
n = cross(J[:,1], J[:,2])
|
||||
n /= norm(n)
|
||||
# sign convention, positive pressure is towards surface
|
||||
p = -element("surface pressure", ip, time)
|
||||
f += w*p*vec(n*N)
|
||||
if has_concentrated_forces
|
||||
update!(element, "concentrated force", time => Any[f])
|
||||
end
|
||||
|
||||
gdofs = get_gdofs(problem, element)
|
||||
add!(assembly.f, gdofs, f)
|
||||
|
||||
end
|
||||
if has_concentrated_forces
|
||||
update!(element, "concentrated force", time => Any[f])
|
||||
end
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
""" Return strain tensor. """
|
||||
|
||||
@@ -4,7 +4,9 @@
|
||||
const Elasticity2DSurfaceElements = Union{Poi1,Seg2,Seg3}
|
||||
const Elasticity2DVolumeElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9}
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_stress}})
|
||||
function assemble!{T}(assembly::Assembly, problem::Problem{Elasticity},
|
||||
elements::Union{Vector{Element}, Vector{Element{T}}},
|
||||
time, ::Type{Val{:plane_stress}})
|
||||
for element in elements
|
||||
gdofs = get_gdofs(problem, element)
|
||||
Km, Kg, f = assemble(problem, element, time, Val{:plane})
|
||||
@@ -14,7 +16,9 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::V
|
||||
end
|
||||
end
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_strain}})
|
||||
function assemble!{T}(assembly::Assembly, problem::Problem{Elasticity},
|
||||
elements::Union{Vector{Element}, Vector{Element{T}}},
|
||||
time, ::Type{Val{:plane_strain}})
|
||||
for element in elements
|
||||
gdofs = get_gdofs(problem, element)
|
||||
Km, Kg, f = assemble(problem, element, time, Val{:plane})
|
||||
@@ -25,7 +29,9 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::V
|
||||
end
|
||||
|
||||
""" Plane elasticity equations (plane stress, plane strain). """
|
||||
function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
|
||||
function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
|
||||
element::Element{El}, time,
|
||||
::Type{Val{:plane}})
|
||||
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
@@ -178,7 +184,9 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
|
||||
return Km, Kg, f
|
||||
end
|
||||
|
||||
function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
|
||||
function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity},
|
||||
element::Element{El},
|
||||
time, ::Type{Val{:plane}})
|
||||
|
||||
props = problem.properties
|
||||
dim = get_unknown_field_dimension(problem)
|
||||
|
||||
Reference in New Issue
Block a user