Skip to content

Commit

Permalink
Manual changes after Runic formatting
Browse files Browse the repository at this point in the history
  • Loading branch information
fredrikekre committed Oct 29, 2024
1 parent 9bd3068 commit 21fec30
Show file tree
Hide file tree
Showing 40 changed files with 1,360 additions and 1,369 deletions.
3 changes: 1 addition & 2 deletions docs/src/literate-gallery/helmholtz.jl
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,7 @@ update!(dbcs, 0.0)
K = allocate_matrix(dh);

function doassemble(
cellvalues::CellValues, facetvalues::FacetValues,
K::SparseMatrixCSC, dh::DofHandler
cellvalues::CellValues, facetvalues::FacetValues, K::SparseMatrixCSC, dh::DofHandler
)
b = 1.0
f = zeros(ndofs(dh))
Expand Down
16 changes: 8 additions & 8 deletions docs/src/literate-gallery/landau.jl
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,7 @@ using Base.Threads
# ### 4th order Landau free energy
function Fl(P::Vec{3, T}, α::Vec{3}) where {T}
P2 = Vec{3, T}((P[1]^2, P[2]^2, P[3]^2))
return (
α[1] * sum(P2) +
α[2] * (P[1]^4 + P[2]^4 + P[3]^4)
) +
return (α[1] * sum(P2) + α[2] * (P[1]^4 + P[2]^4 + P[3]^4)) +
α[3] * ((P2[1] * P2[2] + P2[2] * P2[3]) + P2[1] * P2[3])
end
# ### Ginzburg free energy
Expand Down Expand Up @@ -116,9 +113,10 @@ end

# utility to quickly save a model
function save_landau(path, model, dofs = model.dofs)
return VTKGridFile(path, model.dofhandler) do vtk
VTKGridFile(path, model.dofhandler) do vtk
write_solution(vtk, model.dofhandler, dofs)
end
return
end

# ## Assembly
Expand Down Expand Up @@ -161,19 +159,21 @@ end
# The gradient calculation for each dof
function ∇F!(∇f::Vector{T}, dofvector::Vector{T}, model::LandauModel{T}) where {T}
fill!(∇f, zero(T))
return @assemble! begin
@assemble! begin
ForwardDiff.gradient!(cache.element_gradient, cache.element_potential, eldofs, cache.gradconf)
@inbounds assemble!(∇f, cache.element_indices, cache.element_gradient)
end
return
end

# The Hessian calculation for the whole grid
function ∇²F!(∇²f::SparseMatrixCSC, dofvector::Vector{T}, model::LandauModel{T}) where {T}
assemblers = [start_assemble(∇²f) for t in 1:nthreads()]
return @assemble! begin
@assemble! begin
ForwardDiff.hessian!(cache.element_hessian, cache.element_potential, eldofs, cache.hessconf)
@inbounds assemble!(assemblers[threadid()], cache.element_indices, cache.element_hessian)
end
return
end

# We can also calculate all things in one go!
Expand Down Expand Up @@ -204,7 +204,7 @@ function minimize!(model; kwargs...)
end
function h!(storage, x)
return ∇²F!(storage, x, model)
#apply!(storage, model.boundaryconds)
# apply!(storage, model.boundaryconds)
end
f(x) = F(x, model)

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -262,8 +262,8 @@ function assemble_global!(
## Loop over all cells in the grid
for cell in CellIterator(dh)
global_dofs = celldofs(cell)
global_dofsu = global_dofs[1:nu] # first nu dofs are displacement
global_dofsp = global_dofs[(nu + 1):end] # last np dofs are pressure
global_dofsu = global_dofs[1:nu] # first nu dofs are displacement
global_dofsp = global_dofs[(nu + 1):end] # last np dofs are pressure
@assert size(global_dofs, 1) == nu + np # sanity check
ue = w[global_dofsu] # displacement dofs for the current cell
pe = w[global_dofsp] # pressure dofs for the current cell
Expand Down
24 changes: 12 additions & 12 deletions docs/src/literate-gallery/topology_optimization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,7 @@ function cache_neighborhood(dh, topology)
i = cellid(element)
for j in 1:_nfacets
nbg_cellid = getneighborhood(topology, dh.grid, FacetIndex(i, j))
if (!isempty(nbg_cellid))
if !isempty(nbg_cellid)
nbg[j] = first(nbg_cellid)[1] # assuming only one facet neighbor per cell
else # boundary facet
nbg[j] = first(getneighborhood(topology, dh.grid, FacetIndex(i, opp[j])))[1]
Expand Down Expand Up @@ -246,7 +246,7 @@ function compute_χn1(χn, Δχ, ρ, ηs, χ_min)
λ_upper = maximum(Δχ) + ηs
λ_trial = 0.0

while (abs- ρ_trial) > 1.0e-7)
while abs- ρ_trial) > 1.0e-7
for i in 1:n_el
Δχt = 1 / ηs * (Δχ[i] - λ_trial)
χ_trial[i] = max(χ_min, min(1.0, χn[i] + Δχt))
Expand All @@ -257,9 +257,9 @@ function compute_χn1(χn, Δχ, ρ, ηs, χ_min)
ρ_trial += χ_trial[i] / n_el
end

if (ρ_trial > ρ)
if ρ_trial > ρ
λ_lower = λ_trial
elseif (ρ_trial < ρ)
elseif ρ_trial < ρ
λ_upper = λ_trial
end
λ_trial = 1 / 2 * (λ_upper + λ_lower)
Expand Down Expand Up @@ -304,7 +304,7 @@ function update_density(dh, states, mp, ρ, neighboorhoods, Δh)

χn1 = compute_χn1(χn, Δχ, ρ, mp.η, mp.χ_min)

if (j < n_j)
if j < n_j
χn[:] = χn1[:]
end
end
Expand Down Expand Up @@ -364,7 +364,7 @@ function elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)

symmetrize_lower!(Ke)

return @inbounds for facet in 1:nfacets(getcells(grid, cellid(element)))
@inbounds for facet in 1:nfacets(getcells(grid, cellid(element)))
if (cellid(element), facet) getfacetset(grid, "traction")
reinit!(facetvalues, element, facet)
t = Vec((0.0, -1.0)) # force pointing downwards
Expand All @@ -377,7 +377,7 @@ function elmt!(Ke, re, element, cellvalues, facetvalues, grid, mp, ue, state)
end
end
end

return
end

function symmetrize_lower!(K)
Expand Down Expand Up @@ -479,13 +479,13 @@ function topopt(ra, ρ, n, filename; output = :false)
## calculate compliance
compliance = 1 / 2 * u' * K * u

if (it == 1)
if it == 1
compliance_0 = compliance
end

## check convergence criterium (twice!)
if (abs(compliance - compliance_n) / compliance < tol)
if (conv)
if abs(compliance - compliance_n) / compliance < tol
if conv
println("Converged at iteration number: ", it)
break
else
Expand All @@ -505,7 +505,7 @@ function topopt(ra, ρ, n, filename; output = :false)
compliance_n = compliance

## output during calculation
if (output)
if output
i = @sprintf("%3.3i", it)
filename_it = string(filename, "_", i)

Expand All @@ -516,7 +516,7 @@ function topopt(ra, ρ, n, filename; output = :false)
end

## export converged results
if (!output)
if !output
VTKGridFile(filename, grid) do vtk
write_cell_data(vtk, χ, "density")
end
Expand Down
1 change: 1 addition & 0 deletions docs/src/literate-howto/threaded_assembly.jl
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,7 @@ function create_example_2d_grid()
Ferrite.write_cell_colors(vtk, grid, colors_workstream, "workstream-coloring")
Ferrite.write_cell_colors(vtk, grid, colors_greedy, "greedy-coloring")
end
return
end

create_example_2d_grid()
Expand Down
36 changes: 16 additions & 20 deletions docs/src/literate-tutorials/computational_homogenization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -458,29 +458,25 @@ end
# So we have now already computed all the components, and just need to gather the data in
# a fourth order tensor:

E_dirichlet = SymmetricTensor{4, 2}(
(i, j, k, l) -> begin
if k == l == 1
σ̄.dirichlet[1][i, j] # ∂σ∂ε_**11
elseif k == l == 2
σ̄.dirichlet[2][i, j] # ∂σ∂ε_**22
else
σ̄.dirichlet[3][i, j] # ∂σ∂ε_**12 and ∂σ∂ε_**21
end
E_dirichlet = SymmetricTensor{4, 2}() do i, j, k, l
if k == l == 1
σ̄.dirichlet[1][i, j] # ∂σ∂ε_**11
elseif k == l == 2
σ̄.dirichlet[2][i, j] # ∂σ∂ε_**22
else
σ̄.dirichlet[3][i, j] # ∂σ∂ε_**12 and ∂σ∂ε_**21
end
)
end

E_periodic = SymmetricTensor{4, 2}(
(i, j, k, l) -> begin
if k == l == 1
σ̄.periodic[1][i, j]
elseif k == l == 2
σ̄.periodic[2][i, j]
else
σ̄.periodic[3][i, j]
end
E_periodic = SymmetricTensor{4, 2}() do i, j, k, l
if k == l == 1
σ̄.periodic[1][i, j]
elseif k == l == 2
σ̄.periodic[2][i, j]
else
σ̄.periodic[3][i, j]
end
);
end

# We can check that the result are what we expect, namely that the stiffness with Dirichlet
# boundary conditions is higher than when using periodic boundary conditions, and that
Expand Down
3 changes: 2 additions & 1 deletion docs/src/literate-tutorials/hyperelasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -298,12 +298,13 @@ function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN)
assembler = start_assemble(K, g)

## Loop over all cells in the grid
return @timeit "assemble" for cell in CellIterator(dh)
@timeit "assemble" for cell in CellIterator(dh)
global_dofs = celldofs(cell)
ue = u[global_dofs] # element dofs
@timeit "element assemble" assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN)
assemble!(assembler, global_dofs, ke, ge)
end
return
end;

# Finally, we define a main function which sets up everything and then performs Newton
Expand Down
3 changes: 2 additions & 1 deletion docs/src/literate-tutorials/incompressible_elasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -274,7 +274,8 @@ function solve(ν, interpolation_u, interpolation_p)
σvM = map(x -> (3 / 2 * dev(x) dev(x)), σ) # von Mise effective stress

## Export the solution and the stress
filename = "cook_" * (interpolation_u == Lagrange{RefTriangle, 1}()^2 ? "linear" : "quadratic") *
filename = "cook_" *
(interpolation_u == Lagrange{RefTriangle, 1}()^2 ? "linear" : "quadratic") *
"_linear"

VTKGridFile(filename, grid) do vtk
Expand Down
4 changes: 2 additions & 2 deletions docs/src/literate-tutorials/linear_elasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -358,8 +358,8 @@ colors = [ #hide
"1" => 1, "5" => 1, # purple #hide
"2" => 2, "3" => 2, # red #hide
"4" => 3, # blue #hide
"6" => 4, # green #hide
] #hide
"6" => 4, # green #hide
] #hide
for (key, color) in colors #hide
for i in getcellset(grid, key) #hide
color_data[i] = color #hide
Expand Down
3 changes: 2 additions & 1 deletion docs/src/literate-tutorials/linear_shell.jl
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,8 @@ function fiber_coordsys(Ps::Vector{Vec{3, Float64}})
a = abs.(P)
j = 1
if a[1] > a[3]
a[3] = a[1]; j = 2
a[3] = a[1]
j = 2
end
if a[2] > a[3]
j = 3
Expand Down
27 changes: 12 additions & 15 deletions docs/src/literate-tutorials/plasticity.jl
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticit
σᵗ = material.Dᵉ - state.ϵᵖ) # trial-stress
sᵗ = dev(σᵗ) # deviatoric part of trial-stress
J₂ = 0.5 * sᵗ sᵗ # second invariant of sᵗ
σᵗₑ = sqrt(3.0 * J₂) # effective trial-stress (von Mises stress)
σᵗₑ = sqrt(3.0 * J₂) # effective trial-stress (von Mises stress)
σʸ = material.σ₀ + H * state.k # Previous yield limit

φᵗ = σᵗₑ - σʸ # Trial-value of the yield surface
Expand All @@ -109,7 +109,7 @@ function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticit
return σᵗ, material.Dᵉ, MaterialState(state.ϵᵖ, σᵗ, state.k)
else # plastic loading
h = H + 3G
μ = φᵗ / h # plastic multiplier
μ = φᵗ / h # plastic multiplier

c1 = 1 - 3G * μ / σᵗₑ
s = c1 * sᵗ # updated deviatoric stress
Expand All @@ -129,8 +129,8 @@ function compute_stress_tangent(ϵ::SymmetricTensor{2, 3}, material::J2Plasticit

## Return new state
Δϵᵖ = 3 / 2 * μ / σₑ * s # plastic strain
ϵᵖ = state.ϵᵖ + Δϵᵖ # plastic strain
k = state.k + μ # hardening variable
ϵᵖ = state.ϵᵖ + Δϵᵖ # plastic strain
k = state.k + μ # hardening variable
return σ, D, MaterialState(ϵᵖ, σ, k)
end
end
Expand Down Expand Up @@ -200,10 +200,7 @@ end
#md # Due to symmetry, we only compute the lower half of the tangent
#md # and then symmetrize it.
#md #
function assemble_cell!(
Ke, re, cell, cellvalues, material,
ue, state, state_old
)
function assemble_cell!(Ke, re, cell, cellvalues, material, ue, state, state_old)
n_basefuncs = getnbasefunctions(cellvalues)
reinit!(cellvalues, cell)

Expand All @@ -222,7 +219,8 @@ function assemble_cell!(
end
end
end
return symmetrize_lower!(Ke)
symmetrize_lower!(Ke)
return
end

# Helper function to symmetrize the material tangent
Expand Down Expand Up @@ -257,10 +255,10 @@ end
# Define a function which solves the FE-problem.
function solve()
## Define material parameters
E = 200.0e9 # [Pa]
E = 200.0e9 # [Pa]
H = E / 20 # [Pa]
ν = 0.3 # [-]
σ₀ = 200.0e6 # [Pa]
ν = 0.3 # [-]
σ₀ = 200.0e6 # [Pa]
material = J2Plasticity(E, ν, σ₀, H)

L = 10.0 # beam length [m]
Expand All @@ -285,10 +283,10 @@ function solve()

## Pre-allocate solution vectors, etc.
n_dofs = ndofs(dh) # total number of dofs
u = zeros(n_dofs) # solution vector
u = zeros(n_dofs) # solution vector
Δu = zeros(n_dofs) # displacement correction
r = zeros(n_dofs) # residual
K = allocate_matrix(dh) # tangent stiffness matrix
K = allocate_matrix(dh) # tangent stiffness matrix

## Create material states. One array for each cell, where each element is an array of material-
## states - one for each integration point
Expand All @@ -310,7 +308,6 @@ function solve()

while true
newton_itr += 1

if newton_itr > 8
error("Reached maximum Newton iterations, aborting")
break
Expand Down
3 changes: 2 additions & 1 deletion docs/src/literate-tutorials/porous_media.jl
Original file line number Diff line number Diff line change
Expand Up @@ -359,7 +359,8 @@ function solve(dh, ch, domains; Δt = 0.025, t_total = 1.0)
pvd[t] = vtk
end
end
return vtk_save(pvd)
vtk_save(pvd)
return
end;

# Finally we call the functions to actually run the code
Expand Down
9 changes: 5 additions & 4 deletions docs/src/literate-tutorials/reactive_surface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,8 @@ function setup_initial_conditions!(u₀::Vector, cellvalues::CellValues, dh::Dof
end
end

return u₀ .+= 0.01 * rand(ndofs(dh))
u₀ .+= 0.01 * rand(ndofs(dh))
return
end;

# ### Mesh generation
Expand All @@ -227,7 +228,7 @@ function create_embedded_sphere(refinements)
nodes = tonodes()
elements, _ = toelements(2)
gmsh.finalize()
return grid = Grid(elements, nodes)
return Grid(elements, nodes)
end

# ### Simulation routines
Expand Down Expand Up @@ -312,8 +313,8 @@ function gray_scott_on_sphere(material::GrayScottMaterial, Δt::Real, T::Real, r
## Finally we totate the solution to initialize the next timestep.
uₜ₋₁ .= uₜ
end

return vtk_save(pvd)
vtk_save(pvd)
return
end

## This parametrization gives the spot pattern shown in the gif above.
Expand Down
Loading

0 comments on commit 21fec30

Please sign in to comment.