Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-26 11:51:31 +00:00
Code Issues Packages Projects Releases Wiki Activity

yet another attempt to solve overdetermined system. plane strain formulation.

Browse Source
This commit is contained in:
Jukka Aho
2016-02-08 02:35:39 +02:00
parent 4a78a66aff
commit 6e90492c7b
2 changed files with 65 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elasticity.jl
+22 -6
View File
@@ -26,12 +26,18 @@ function get_unknown_field_name(::Type{Elasticity})
end
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real)
return assemble!(assembly, problem, element, time, Val{problem.properties.formulation})
f = problem.properties.formulation
if f == :continuum
return assemble!(assembly, problem, element, time, Val{:continuum})
elseif (f == :plane_stress) || (f == :plane_strain)
return assemble!(assembly, problem, element, time, Val{:plane})
end
end
""" Elasticity equations, plane stress formulation. """
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:plane_stress}})
function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:plane}})
props = problem.properties
gdofs = get_gdofs(problem, element)
ndim, nnodes = size(element)
B = zeros(3, 2*nnodes)
@@ -42,10 +48,20 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El
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]
# Zienkiewicz, p. 91
if props.formulation == :plane_stress
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]
elseif props.formulation == :plane_strain
C = 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
dN = element(ip, time, Val{:grad})
fill!(B, 0.0)
for i=1:size(dN, 2)
src/solver_utils.jl
+43 -3
View File
@@ -91,7 +91,15 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
function is_spc(dofs::Vector{Int})
return map(is_spc, dofs)
end
function has_anything(dof::Int)
countnz(C1[dof,:]) != 0 && return true
countnz(C2[dof,:]) != 0 && return true
countnz(D[dof,:]) != 0 && return true
countnz(g[dof,:]) != 0 && return true
return false
end
""" Algorithm 1. Calculate rank of overdetermined system and do LSQ if
rank(C) equals to number of unique dofs.
"""
@@ -147,7 +155,39 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
end
end
actions = [action1, action3]
function action4(node_id, dofs)
""" If symmetry line, one possibility is to apply both conditions and
eliminate lagrange multiplier. """
dofs_ = intersect(dofs, all_dofs)
length(dofs_) != 1 && return dofs_, false
related_dofs = get_related_dofs(dofs_)
for j in related_dofs
has_anything(j) && continue
# copy one constaint to this dof
C1[j,:] = C1[dofs_,:]
C2[j,:] = C2[dofs_,:]
D[j,:] = D[dofs_,:]
g[j,:] = g[dofs_,:]
# make room for new constraint
C1[dofs_,:] = 0
C2[dofs_,:] = 0
D[dofs_,:] = 0
g[dofs_,:] = 0
dofs_ = [dofs_; j]
break
end
for j in related_dofs
has_anything(j) && continue
# set lagrange multiplier to 1
D[j,dofs_[1]] = 1.0
C1[j,dofs_[1]] = 1.0
dofs_ = [dofs_; j]
break
end
return dofs_, true
end
actions = [action1, action2]
function show_lambda_coefficients(dofs, C1)
for dof in dofs
@@ -232,7 +272,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
info("fixed: new setting is")
show_rows_in_constraint_matrix(dofs, C2, D; show_status=false)
show_rows_in_constraint_matrix(dofs, C2_, D_; show_status=false)
show_related_equations(dofs, C2, C2_, D, D_)
#show_related_equations(dofs, C2, C2_, D, D_)
info()
continue
end