Finite Element Analysis with Gridap.jl; how to define external forces?

Question

I'm following this tutorial in order to try and do an FEA of a model.msh that I have to see how it would deform given different external forces in different places.

There they define the weak form as

a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ
l(v) = 0

and they state "The linear form is simply l(v) = 0 since there are not external forces in this example."

As mentioned, I would like to analyse the different deformation that different external forces would cause on my model, but I can't seem to find anywhere an example of this. Could someone help me on defining this linear form for external forces different than 0?

Thanks.

Answer 1

Maybe this helps you out. It was written in a hurry, so please do not mind if you encounter spelling mistakes or other beauty issues :)

# define where the output shall go
output_Path ="Output/3_Gridap/1_Lin_FEA/FE_8"
mkpath(output_Path)
output_Name ="pde_6"

using Gridap

# please load the model that is shown in: https://gridap.github.io/Tutorials/dev/pages/t001_poisson/
model = DiscreteModelFromFile("Code/Meshes/Data/possion.json")

# just in case you want to see the model using paraview
writevtk(model,"$(output_Path)/md")

order = 1
reffe = ReferenceFE(lagrangian,VectorValue{3,Float64},order)

V0 = TestFESpace(model,reffe;
  conformity=:H1,

  # to see which elements belongs to "bottom" open the model which is saved through "writevtk(model,"$(output_Path)/md")"
  dirichlet_tags=["bottom"],
  
  # activate/deactivate the boundary conditions
  dirichlet_masks=[
        (true, true, true),       # clamp the bottom
        ])

# define displacement
clamping(x) = VectorValue(0.0,0.0,0.0)

U = TrialFESpace(V0,[clamping])
const E = 7e+7
const ν = 0.33
const λ = (E*ν)/((1+ν)*(1-2*ν))
const μ = E/(2*(1+ν))
σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε

degree = 2*order
Ω = Triangulation(model)
dΩ = Measure(Ω,degree)

# Neumann boundary conditions
# here we define the surface on which we want an external force
Γ_Tr = BoundaryTriangulation(model,tags=["triangle"])
dΓ_Tr = Measure(Γ_Tr,degree)

# a force shall be applied on the y-direction
f_Tr(x) = VectorValue(0.0, 1e+6, 0.0)

# mass forces due to gravity, the value is set quite high, such that an impact can be seen
mass_Forces(x) = VectorValue(0.0, -1e+7, 0.0)

# Weak form
a(u,v) = ∫( ε(v) ⊙ (σ∘ε(u)) )*dΩ

l(v) = ∫( v ⋅ mass_Forces )* dΩ + ∫( v ⋅ f_Tr )* dΓ_Tr 

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

writevtk(Ω,"$(output_Path)/$(output_Name)",
      
cellfields=[
            "uh" => uh,
            "epsi" => ε(uh),
            "sigma" => σ∘ε(uh)])