 # Assembly prescribed Stress using JuAFEM

Checking the example incompressible_elasticity, it is solving the equation Ku = f based on the following assemble code:

``````function doassemble(cellvalues_u::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim},
facevalues_u::FaceVectorValues{dim}, K::SparseMatrixCSC, grid::Grid,
dh::DofHandler, mp::LinearElasticity) where {dim}

f = zeros(ndofs(dh))
assembler = start_assemble(K, f)
nu = getnbasefunctions(cellvalues_u)
np = getnbasefunctions(cellvalues_p)

fe = PseudoBlockArray(zeros(nu + np), [nu, np]) # local force vector
ke = PseudoBlockArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix

# traction vector
t = Vec{2}((0.0, 1/16))
# cache ɛdev outside the element routine to avoid some unnecessary allocations
ɛdev = [zero(SymmetricTensor{2, dim}) for i in 1:getnbasefunctions(cellvalues_u)]

for cell in CellIterator(dh)
fill!(ke, 0)
fill!(fe, 0)
assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, grid, mp, ɛdev, t)
assemble!(assembler, celldofs(cell), fe, ke)
end

return K, f
end;
``````

And after applying boundary conditions by `apply!(K, f, dbc)` one can solve the system.

My question is, what is the best way to solve the system Ku = f -g, where g is a vector describing predescribed stress. My first was, was to create a vector g_e and subtract it from the vector f_e such that f = Assemble(f_e -g_e). The problem is, that the command `apply!(K, f, dbc)` would apply boundary conditions on g too. And I am not sure if this is the correct way.

Heyho,

IIRC neumann and dirichlet boundary do not intersect. Therefore, I don’t think that `apply!(K,f,dbc)` does anything wrong for you, since g_e should be 0 anyways at the dirichlet boundary. You can test this with: https://github.com/KristofferC/JuAFEM.jl/pull/304

1 Like