Assembly prescribed Stress using JuAFEM

General Usage
1

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.

2

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