Compute Compliance along with multiple force using Gridap

Code

using Gridap
using GridapGmsh
using LinearAlgebra
using Random
using Gridap, Gridap.Geometry, Gridap.TensorValues


model = DiscreteModelFromFile("./model.json")

## 
# 
# 
order = 1
degree = 2*order

reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
V = TestFESpace(model,reffe;
  conformity=:H1,
  dirichlet_tags=["dirichlet"],
  dirichlet_masks=[(true,true,true)]
)
U = TrialFESpace(V)
free_ids = Gridap.FESpaces.get_free_dof_ids(V)
dir_ids  = Gridap.FESpaces.get_dirichlet_dof_ids(V)

println("free dofs = ", length(free_ids))
println("fixed dofs = ", length(dir_ids))
println("total dofs = ", length(free_ids) + length(dir_ids))

Ω  = Triangulation(model)
dΩ = Measure(Ω, degree)


# 
# 
# 
E  = 210e9
ν = 0.3
λ = (E*ν)/((1+ν)*(1-2ν))
μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε


ε(u) = symmetric_gradient(u)
σ_point(εp) = 2μ*εp + λ*tr(εp)*I

const d = 3
const I = one(Gridap.TensorValues.SymTensorValue{d,Float64})
a(u,v) = ∫( (2μ*ε(u) + λ*tr(ε(u))*I) ⊙ ε(v) )dΩ
# a(u,v) = ∫( ε(v) ⊙ (σ_point ∘ ε(u)))dΩ
# a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)))dΩ

# 
# Assemble Stiffness Matrix
# 
A = assemble_matrix(a, U, V)
labels = Gridap.Geometry.get_face_labeling(model)
Gridap.Geometry.add_tag_from_tags_complementary!(labels, "not_dirichlet", ["dirichlet"])
Γ = Gridap.Geometry.BoundaryTriangulation(model, tags="not_dirichlet")
# Γ = Gridap.Geometry.BoundaryTriangulation(model, tags=force_tags)
dΓ = Measure(Γ, degree)

N = 10
nc = Gridap.Geometry.num_cells(Γ)
if nc == 0
    error("There is no facets on the Γ.")
end
@assert nc > 0 "There is no facet on the surface Γ"

k = min(N, nc)
rand_ids = k == 0 ? Int[] : randperm(nc)[1:k]

ndofs = Gridap.FESpaces.num_free_dofs(V)
B = zeros(ndofs, k)

nΓ  = Gridap.CellData.get_normal_vector(Γ)
t0  = 1e6
t_all = t0 * nΓ

for (j, cid) in enumerate(rand_ids)
    mask_vals = [i==cid ? 1.0 : 0.0 for i in 1:nc]
    m = Gridap.CellData.CellField(mask_vals, Γ)
    t_cell = m * t_all
    l(v) = ∫( t_cell ⋅ v )dΓ
    b = assemble_vector(l, V)
    B[:, j] .= b
end


A_lu = lu(A) # sparse lu from UMFPACK, not CHOLMOD
Ucols = A_lu \ B #can be solved simultaniously
C = vec(sum(B .* Ucols; dims=1))  # list of compliance

println("num faces = $nc")
println("compliance list: $C")
println("compliance stats: min=$(minimum(C)), max=$(maximum(C))")