Gridap.jl's Linear Elasticity Problem: Access to Global and Local Stiffness matrixes

General Usage
1

Hello everyone!
I solved the Linear elasticity problem for a cantilever beam.

image
image506×312 57.6 KB

Now I want to compute the compliance of structures as follows:

c = \mathbf{U}^T \mathbf{K} \mathbf{U} = \sum_{e=1}^{N} \mathbf{u}_e^T \mathbf{k}_e\, \mathbf{u}_e.

To do this, I need access to the global stiffness matrix and the degree of freedom of each element (dof).
Could you please help me 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
2

I just found the solutions. Here are three methods:

# 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)
3

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