In [15]:
using JuAFEM
using Tensors
using TimerOutputs
using UnicodePlots
const to = TimerOutput();
In [16]:
hex = false
geoshape = hex ? Hexahedron : Tetrahedron
refshape = hex ? RefCube : RefTetrahedron
order = hex? 2 : 1;
In [17]:
const dim = 3
corner1 = Vec{dim}((0.0, 0.0, 0.0))
corner2 = Vec{dim}((10.0, 1.0, 1.0))
grid = generate_grid(geoshape, (60, 6, 6), corner1, corner2);
# Extract the left boundary
addnodeset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0);
In [18]:
# Interpolations and values
interpolation_space = Lagrange{dim, refshape, 1}()
quadrature_rule = QuadratureRule{dim, refshape}(order)
cellvalues = CellVectorValues(quadrature_rule, interpolation_space);
facevalues = FaceVectorValues(QuadratureRule{dim-1, refshape}(order), interpolation_space);
In [19]:
# DofHandler
dh = DofHandler(grid)
push!(dh, :u, dim) # Add a displacement field
close!(dh)
Out[19]:
In [20]:
@time K = create_symmetric_sparsity_pattern(dh); # assemble only upper half since it is symmetric
fill!(K.data.nzval, 1.0);
spy(K.data)
Out[20]:
In [21]:
# Boundaryconditions
dbc = DirichletBoundaryConditions(dh)
# Add a homogenoush boundary condition on the "clamped" edge
add!(dbc, :u, getnodeset(grid, "clamped"), (x,t) -> [0.0, 0.0, 0.0], collect(1:dim))
close!(dbc)
t = 0.0
update!(dbc, t)
In [22]:
# Create the stiffness tensor
E = 200e9
ν = 0.3
λ = E*ν / ((1+ν) * (1 - 2ν))
μ = E / (2(1+ν))
δ(i,j) = i == j ? 1.0 : 0.0
g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
C = SymmetricTensor{4, dim}(g);
In [23]:
function doassemble{dim}(cellvalues::CellVectorValues{dim}, facevalues::FaceVectorValues{dim},
K::Symmetric, grid::Grid, dh::DofHandler, C::SymmetricTensor{4, dim})
f = zeros(ndofs(dh))
assembler = start_assemble(K, f)
n_basefuncs = getnbasefunctions(cellvalues)
fe = zeros(n_basefuncs) # Local force vector
Ke = Symmetric(zeros(n_basefuncs, n_basefuncs), :U) # Local stiffness mastrix
t = Vec{3}((0.0, 1e8, 0.0)) # Traction vector
b = Vec{3}((0.0, 0.0, 0.0)) # Body force
ɛ = [zero(SymmetricTensor{2, dim}) for i in 1:n_basefuncs]
@inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
@timeit to "assem" begin
fill!(Ke.data, 0)
fill!(fe, 0)
reinit!(cellvalues, cell)
for q_point in 1:getnquadpoints(cellvalues)
for i in 1:n_basefuncs
ɛ[i] = symmetric(shape_gradient(cellvalues, q_point, i))
end
dΩ = getdetJdV(cellvalues, q_point)
for i in 1:n_basefuncs
δu = shape_value(cellvalues, q_point, i)
fe[i] += (δu ⋅ b) * dΩ
ɛC = ɛ[i] ⊡ C
for j in i:n_basefuncs # assemble only upper half
Ke.data[i, j] += (ɛC ⊡ ɛ[j]) * dΩ # can only assign to parent of the Symmetric wrapper
end
end
end
for face in 1:nfaces(cell)
if onboundary(cell, face) && (cellcount, face) ∈ getfaceset(grid, "right")
reinit!(facevalues, cell, face)
for q_point in 1:getnquadpoints(facevalues)
dΓ = getdetJdV(facevalues, q_point)
for i in 1:n_basefuncs
δu = shape_value(facevalues, q_point, i)
fe[i] += (δu ⋅ t) * dΓ
end
end
end
end
global_dofs = celldofs(cell)
assemble!(assembler, global_dofs, fe, Ke)
end # timer
end
return K, f
end;
In [24]:
reset_timer!(to)
K, f = doassemble(cellvalues, facevalues, K, grid, dh, C);
print_timer(to; linechars = :ascii)
In [25]:
# Modify K and f such that K \ f gives correct boundary conditions
@time apply!(K, f, dbc)
In [26]:
@time u = cholfact(K) \ f;
In [27]:
# Save file
vtkfile = vtk_grid("cantilever", dh, u)
vtk_save(vtkfile)
Out[27]:
In [28]:
Base.Test.@test maximum(u) ≈ 1.919600482922295
println("Cantilever successful")
In [ ]: