Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-17 01:02:13 +00:00
Code Issues Packages Projects Releases Wiki Activity

fixed weighted gap in 3d problems + some small fixes. nt formulation

Browse Source
This commit is contained in:
Jukka Aho
2016-02-06 15:29:35 +02:00
parent 070e03eb05
commit 5c5adff939
3 changed files with 40 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elasticity.jl
+2 -2
View File
@@ -143,10 +143,10 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
add!(assembly.f, gdofs, vec(L))
end
for dim in 1:problem.dim
for dim in 1:get_unknown_field_dimension(problem)
if haskey(element, "displacement traction force $dim")
T = element("displacement traction force $dim", ip, time)
ldofs = gdofs[dim:problem.dim:end]
ldofs = gdofs[dim:unknown_field_dimension(problem):end]
JT = transpose(J)
L = w*T*N*norm(cross(JT[:,1], JT[:,2]))
add!(assembly.f, ldofs, vec(L))
src/integrate.jl
+1 -1
View File
@@ -136,7 +136,7 @@ end
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) for i=1:2, j=1:2, k=1:2])
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})
src/mortar.jl
+37 -7
View File
@@ -921,14 +921,44 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mor
Sm = ip.weight*N1'*N1*wC
Mm = ip.weight*N1'*N2*wC
for k=1:field_dim
sd = slave_dofs[k:field_dim:end]
md = master_dofs[k:field_dim:end]
add!(assembly.C1, sd, sd, Sm)
add!(assembly.C1, sd, md, -Mm)
add!(assembly.C2, sd, sd, Sm)
add!(assembly.C2, sd, md, -Mm)
@assert length(slave_dofs) == length(master_dofs)
S3 = zeros(length(slave_dofs), length(slave_dofs))
M3 = zeros(length(master_dofs), length(master_dofs))
Q = slave_element("normal-tangential coordinates", xi_slave, time)
Z = zeros(3, 3)
Q3 = [Q Z Z Z; Z Q Z Z; Z Z Q Z; Z Z Z Q]
#info("N1 = $N1")
#info("N2 = $N2")
#info("size S3 = $(size(S3))")
#info("size M3 = $(size(M3))")
#info("size Sm = $(size(Sm))")
#info("size Mm = $(size(Mm))")
for k=1:field_dim
S3[k:field_dim:end,k:field_dim:end] += Sm
M3[k:field_dim:end,k:field_dim:end] += Mm
end
add!(assembly.C1, slave_dofs, slave_dofs, S3)
add!(assembly.C1, slave_dofs, master_dofs, -M3)
S3 = Q3'*S3
M3 = Q3'*M3
add!(assembly.C2, slave_dofs, slave_dofs, S3)
add!(assembly.C2, slave_dofs, master_dofs, -M3)
X1 = slave_element("geometry", xi_slave, time)
X2 = master_element("geometry", xi_master, time)
g = norm(X2-X1)
#T = transpose(get_jacobian(slave_element, xi_slave, time))
#W = ip.weight*
#info("hard gap = $g, wS = $wS, wM = $wM, wC = $wC")
gh = ip.weight*N1*g*wC
add!(assembly.g, slave_dofs[1:field_dim:end], gh)
#for k=1:field_dim
# sd = slave_dofs[k:field_dim:end]
# md = master_dofs[k:field_dim:end]
# add!(assembly.C1, sd, sd, Sm)
# add!(assembly.C1, sd, md, -Mm)
# add!(assembly.C2, sd, sd, Sm)
# add!(assembly.C2, sd, md, -Mm)
#end
end
# info("breaking on first")
# break