Access to Global and Local Stiffness matrixes

[sy-nguyen-van]

Hello everyone!
I solved the Linear elasticity problem for a cantilever beam. Now I want to compute the compliance of structures as follows:

c = Uᵀ K U = Σₑ uₑᵀ kₑ uₑ

To do this, I need access to the global stiffness matrix and the degree of freedom of each element (dof).
Could you please help with how to do this? Thank you so much.
FYI, here is my code:

using Gridap
L = 120 # Length
W = 24 # Width and Height
domain = (0, L, 0, W, 0, W) # sizes along X, Y, Z
partition = (120, 24, 24) # mesh sizes
model = CartesianDiscreteModel(domain, partition)
# Tag the left face (x = 0); right face (x = L)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"left", [1,3,5,7,13,15,17,19,25]) # left face: corners (1,3,5,7); edges (13,15,17,19); others (25) 
add_tag_from_tags!(labels, "right", [2,4,6,8,14,16,18,20,26]) # Right face: corners (2,4,6,8); edges (14,16,18,20); others (26) 
# Integrating on the domain Ω
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)
# neumann boundaries
Γ = BoundaryTriangulation(model,tags = "right")
dΓ = Measure(Γ,degree)
# Vector-valued FE space
order = 1
reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
Vₕ = TestFESpace(Ω, reffe; conformity=:H1, dirichlet_tags = "left")
Uₕ = TrialFESpace(Vₕ)
# Constitutive law for LINEAR ELASTICITY 
const E = 119e3 # Young's modulus
const ν = 0.3   # Poisson's ratio
const λ = (E*ν)/((1+ν)*(1-2*ν)) # Lamé parameters
const μ = E/(2*(1+ν)) # Lamé parameters
ε₀(u) = 0.5 * ( ∇(u) + transpose(∇(u))) # Strain; In Gridap this can be automatically defined in Gridap.ε
σ(ε₀) = λ * tr(ε₀) * one(ε₀) + 2μ * ε₀ # Stress
# The weak form
F(x) = VectorValue(0.0, -1, 0.0) # I want to apply a unit load along Y axis to each node on the right face
a(u,v) = ∫((σ∘ε₀(u)) ⊙ ε₀(v)) *dΩ # Left-hand size; (∘) Composite functions
l(v) = ∫(F ⋅ v) * dΓ # Right-hand size
# Solution of the FE problem
op = AffineFEOperator(a,l,Uₕ,Vₕ)
uh = solve(op)
# Export the displacement, strain, stress
writevtk(Ω,"Cantilever_Beam_3D_uh",cellfields=["uh"=>uh,"epsi"=>ε₀(uh),"sigma"=>σ∘ε₀(uh)])
# ========ACCESS TO GLOBAL STIFFNESS MATRIXES =============================
# Compute the compliance 
num_ele = num_cells(Ω) # No of elements
c = 0 # compliance
K = ... # Global stiffness matrix
for e=1:num_ele
    dof_e = .... # degree of freedom of element "e"
    ue = uh(dof_e) # the element displacement vector
    ke = K(dof_e,dof_e) # the element stiffness matrix
    c = c+ transpose(ue)*ke*ue # update compliance
end

[sy-nguyen-van]

Hello everyone, solutions are found. Here is my code:

# Method 1
C1 = sum(∫((σ ∘ ε₀(uh)) ⊙ ε₀(uh)) * dΩ) # Compliance
println("Compliance 1: ", C1)
# Method 2
C2 = sum(∫(F ⋅ uh) * dΓ) # Compliance
println("Compliance 2: ", C2)
# Method 3
u = get_free_dof_values(uh) # displacement vector
K = get_matrix(op) # stiffness matrix
C3 = dot(u, K * u)
println("Compliance 3: ", C3)

[sy-nguyen-van]

For your information, I still want to manually calculate the compliance using the element stiffness matrix and the element displacement vector. However, my results aren’t correct. I suspect the problem might be related to the number of degrees of freedom. Any hints or suggestions to help resolve this would be greatly appreciated. Thank you!

# Method 4 using for loop to compute sum(ue^T*ke*ue)
U = get_free_dof_values(uh) # displacement vector
K = get_matrix(op) # stiffness matrix
global_dof = get_cell_dof_ids(Vₕ) # Degree of freedom
numele = num_cells(Ω) # Number of elements
C4 = 0.0 # Comliance
for e=1:numele
    index_e = findall(global_dof[e] .>0) # Somehow in Gridap.jl dof of dirichlet_tags is negative
    dof_e = global_dof[e][index_e] # none negative dof
    ue = U[dof_e]    # element displacement
    ke = K[dof_e,dof_e] # element stiffness matrix
    C4 = C4 + transpose(ue)*ke*ue
end
println("Compliance 4: ", C4)
# ---
Compliance 1: 0.30855010386860376
Compliance 2: 0.30855010386853077
Compliance 3: 0.30855010386865317
Compliance 4: 5808.737629988715

[JordiManyer]

Your loop is wrong in two ways: First, this

K = get_matrix(op)

gives you the global assembled stiffness matrix. The restriction of the assembled matrix to the dofs of a cell is NOT the local stiffness matrix. The local stiffness matrices are obtained BEFORE ASSEMBLY, i.e

using Gridap.Arrays
u = get_trial_fe_basis(V)
v = get_fe_basis(V)
K_e = get_array(∫((σ ∘ ε₀(u)) ⊙ ε₀(v)) * dΩ)

Second, you are forgetting about the contributions from the Dirichlet dofs, which you can't. The complete cell dofs you can get as

using Gridap.FESpaces
cell_dofs = get_cell_dof_values(uh)

That said, just use an integral like you did in your first post.

[sy-nguyen-van]

Thank you so much for your information. Great! It works now!

using Gridap.Arrays, Gridap.FESpaces
cell_dofs = get_cell_dof_values(uh)
u = get_trial_fe_basis(V)
v = get_fe_basis(V)
K_e = get_array(∫((σ ∘ ε₀(u)) ⊙ ε₀(v)) * dΩ)
numele = num_cells(Ω) # Number of elements
C4 = 0.0 # Comliance
for e=1:numele
    ue = cell_dofs[e]    # element displacement
    ke = K_e[e] # element stiffness matrix
    C4 = C4 + transpose(ue)*ke*ue
end
println("Compliance 4: ", C4)

Compliance 1: 0.30855010386860376
Compliance 2: 0.30855010386853077
Compliance 3: 0.30855010386865317
Compliance 4: 0.30855010386887693

[sy-nguyen-van]

Dear @JordiManyer,

1. I just wonder in Gridap.jl, are there any low-API functions to compute "the symmetric gradient of the shape functions (B_e)" and "the elasticity tensor (C_e)" as shown in the attached figure ? The reason for accessing C_e and B_e is that I want to compute the stress for each element manually.

# ==Low API function to compute Be and Ce
U_cells = get_cell_dof_values(uh) # displacement vector
numele = num_cells(Ω) # Number of elements
for e=1:numele
    Be = ...  # gradient of the shape functions
    Ce = ... # elasticity tensor
    ue = U_cells[e]    # element displacement
    sigma_e = Ce*Be*ue
end

2. In Gridap.jl, we can get stress by using sigma_cells=σ∘ε₀(uh), but do you know how to get the numeric values of sigma_cells ?

Thanks for your help.

[JordiManyer]

You can get the quadrature points by calling 'get_cell_points(quad)' on the CellQuadrature that the Measure holds.

If you want to learn how the code works internally, I would encourage you to look at the low-level API tutorial on the poisson problem, where all of these things are detailed. You will soon realise what you are asking is a slight variation of the previous answer.
Let me know if you still have issues finding the answer on your own! 😉

[sy-nguyen-van]

Thanks a lot for your hint! I will look at this.
Best regards,