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.