Skip to content

Commit e7c1a79

Browse files
cmhamelcmhamel
authored
Merge pull request #100 from Cthonios/physics/block-by-block
Physics/block by block
2 parents 1434849 + fba56d5 commit e7c1a79

17 files changed

+375
-90
lines changed

‎examples/mechanics/Cube.jl‎

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
using Exodus
2+
using FiniteElementContainers
3+
using Tensors
4+
5+
fixed(_, _) = 0.
6+
displace(_, t) = 1.0 * t
7+
8+
struct Mechanics{Form} <: AbstractPhysics{2, 0, 0}
9+
formulation::Form
10+
end
11+
12+
function strain_energy(∇u)
13+
E = 1.e9
14+
ν = 0.25
15+
K = E / (3. * (1. - 2. * ν))
16+
G = E / (2. * (1. + ν))
17+
18+
ε = symmetric(∇u)
19+
ψ = 0.5 * K * tr(ε)^2 + G * dcontract(dev(ε), dev(ε))
20+
end
21+
22+
function FiniteElementContainers.residual(physics::Mechanics, cell, u_el, args...)
23+
(; X_q, N, ∇N_X, JxW) = cell
24+
25+
# kinematics
26+
∇u_q = u_el * ∇N_X
27+
∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
28+
# constitutive
29+
P_q = gradient(strain_energy, ∇u_q)
30+
# turn into voigt notation
31+
P_q = extract_stress(physics.formulation, P_q)
32+
G_q = discrete_gradient(physics.formulation, ∇N_X)
33+
f_q = G_q * P_q
34+
return JxW * f_q[:]
35+
# return f_q
36+
end
37+
38+
function FiniteElementContainers.stiffness(physics::Mechanics, cell, u_el, args...)
39+
(; X_q, N, ∇N_X, JxW) = cell
40+
41+
# kinematics
42+
∇u_q = u_el * ∇N_X
43+
∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
44+
# constitutive
45+
K_q = hessian(z -> strain_energy(z), ∇u_q)
46+
# turn into voigt notation
47+
K_q = extract_stiffness(physics.formulation, K_q)
48+
G_q = discrete_gradient(physics.formulation, ∇N_X)
49+
return JxW * G_q * K_q * G_q'
50+
end
51+
52+
function mechanics_test()
53+
mesh = UnstructuredMesh("cube.g")
54+
V = FunctionSpace(mesh, H1Field, Lagrange)
55+
physics = (;
56+
cube=Mechanics(ThreeDimensional())
57+
)
58+
59+
u = VectorFunction(V, :displ)
60+
asm = SparseMatrixAssembler(H1Field, u)
61+
62+
dbcs = DirichletBC[
63+
DirichletBC("displ_x", "ssx-", fixed),
64+
DirichletBC("displ_y", "ssy-", fixed),
65+
DirichletBC("displ_z", "ssz-", fixed),
66+
DirichletBC("displ_z", "ssz+", displace),
67+
]
68+
69+
# pre-setup some scratch arrays
70+
times = TimeStepper(0., 1., 10)
71+
p = create_parameters(asm, physics; dirichlet_bcs=dbcs, times=times)
72+
73+
solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
74+
integrator = QuasiStaticIntegrator(solver)
75+
pp = PostProcessor(mesh, "cube.e", u)
76+
77+
write_times(pp, 1, p.times.time_current[1])
78+
write_field(pp, 1, p.h1_field)
79+
80+
for n in 1:10
81+
evolve!(integrator, p)
82+
write_times(pp, n, p.times.time_current[1])
83+
write_field(pp, n, p.h1_field)
84+
end
85+
close(pp)
86+
end
87+
88+
mechanics_test()

‎examples/mechanics/cube-line-search.yaml‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
type: single
2+
input mesh file: cube.g
3+
output mesh file: cube.e
4+
Exodus output interval: 1
5+
CSV output interval: 0
6+
model:
7+
type: solid mechanics
8+
material:
9+
blocks:
10+
cube: elastic
11+
elastic:
12+
model: linear elastic
13+
elastic modulus: 1.0e+09
14+
Poisson's ratio: 0.25
15+
density: 1000.0
16+
time integrator:
17+
type: quasi static
18+
initial time: 0.0
19+
final time: 1.0
20+
time step: 0.1
21+
boundary conditions:
22+
Dirichlet:
23+
- node set: nsx-
24+
component: x
25+
function: "0.0"
26+
- node set: nsy-
27+
component: y
28+
function: "0.0"
29+
- node set: nsz-
30+
component: z
31+
function: "0.0"
32+
- node set: nsz+
33+
component: z
34+
function: "1.0 * t"
35+
solver:
36+
type: Hessian minimizer
37+
step: full Newton
38+
use line search: true
39+
line search backtrack factor: 0.5
40+
line search decrease factor: 1.0e-04
41+
line search maximum iterations: 16
42+
minimum iterations: 1
43+
maximum iterations: 16
44+
relative tolerance: 1.0e-12
45+
absolute tolerance: 1.0e-08

‎examples/mechanics/cube-steepest-descent.yaml‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
type: single
2+
input mesh file: cube.g
3+
output mesh file: cube.e
4+
Exodus output interval: 1
5+
CSV output interval: 0
6+
model:
7+
type: solid mechanics
8+
material:
9+
blocks:
10+
cube: elastic
11+
elastic:
12+
model: linear elastic
13+
elastic modulus: 1.0e+09
14+
Poisson's ratio: 0.25
15+
density: 1000.0
16+
time integrator:
17+
type: quasi static
18+
initial time: 0.0
19+
final time: 1.0
20+
time step: 0.1
21+
boundary conditions:
22+
Dirichlet:
23+
- node set: nsx-
24+
component: x
25+
function: "0.0"
26+
- node set: nsy-
27+
component: y
28+
function: "0.0"
29+
- node set: nsz-
30+
component: z
31+
function: "0.0"
32+
- node set: nsz+
33+
component: z
34+
function: "1.0 * t"
35+
solver:
36+
type: steepest descent
37+
step: steepest descent
38+
step length: 1.0
39+
minimum iterations: 1
40+
maximum iterations: 1024
41+
relative tolerance: 1.0e-15
42+
absolute tolerance: 1.0e-08

‎examples/mechanics/cube.g‎

10.2 KB
Binary file not shown.

‎examples/mechanics/cube.yaml‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
type: single
2+
input mesh file: cube.g
3+
output mesh file: cube.e
4+
Exodus output interval: 1
5+
CSV output interval: 0
6+
model:
7+
type: solid mechanics
8+
material:
9+
blocks:
10+
cube: elastic
11+
elastic:
12+
model: linear elastic
13+
elastic modulus: 1.0e+09
14+
Poisson's ratio: 0.25
15+
density: 1000.0
16+
time integrator:
17+
type: quasi static
18+
initial time: 0.0
19+
final time: 1.0
20+
time step: 0.1
21+
boundary conditions:
22+
Dirichlet:
23+
- node set: nsx-
24+
component: x
25+
function: "0.0"
26+
- node set: nsy-
27+
component: y
28+
function: "0.0"
29+
- node set: nsz-
30+
component: z
31+
function: "0.0"
32+
- node set: nsz+
33+
component: z
34+
function: "1.0 * t"
35+
solver:
36+
type: Hessian minimizer
37+
step: full Newton
38+
minimum iterations: 1
39+
maximum iterations: 16
40+
relative tolerance: 1.0e-12
41+
absolute tolerance: 1.0e-08

‎src/DofManagers.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -219,6 +219,8 @@ Base.length(dof::DofManager) = length(dof.H1_bc_dofs) + length(dof.H1_unknown_do
219219

220220
KA.get_backend(dof::DofManager) = KA.get_backend(dof.H1_unknown_dofs)
221221

222+
# is create_bcs even really needed anymore?
223+
222224
"""
223225
$(TYPEDSIGNATURES)
224226
"""

‎src/FiniteElementContainers.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -143,6 +143,7 @@ include("physics/Physics.jl")
143143
include("PostProcessors.jl")
144144

145145
#
146+
include("TimeSteppers.jl")
146147
include("Parameters.jl")
147148
include("solvers/Solvers.jl")
148149
include("integrators/Integrators.jl")

‎src/Parameters.jl‎

Lines changed: 12 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -1,29 +1,3 @@
1-
# move to a seperate file
2-
abstract type AbstractTimeStepper{T} end
3-
current_time(t::AbstractTimeStepper) = sum(t.time_current)
4-
time_step(t::AbstractTimeStepper) = sum(t.Δt)
5-
6-
struct TimeStepper{T} <: AbstractTimeStepper{T}
7-
time_start::T
8-
time_end::T
9-
time_current::T
10-
Δt::T
11-
end
12-
13-
function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
14-
time_start = zeros(1)
15-
time_end = zeros(1)
16-
time_current = zeros(1)
17-
Δt = zeros(1)
18-
Δt = zeros(1)
19-
fill!(time_start, time_start_in)
20-
fill!(time_end, time_end_in)
21-
fill!(time_current, time_start_in)
22-
fill!(Δt, (time_end_in - time_start_in) / n_steps)
23-
return TimeStepper(time_start, time_end, time_current, Δt)
24-
# return TimeStepper(time_start_in, time_end_in, )
25-
end
26-
271
abstract type AbstractParameters end
282

293
struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
@@ -42,26 +16,13 @@ struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
4216
h1_field::H1
4317
end
4418

45-
# TODO only works for H1Fields currently most likely
46-
# function Parameters(dof::DofManager, physics::AbstractPhysics)
47-
# n_props = num_properties(physics)
48-
# n_state = num_states(physics)
49-
50-
# # TODO
51-
# fspace = dof.H1_vars[1].fspace
52-
53-
# for (name, ref_fe) in pairs(fspace.ref_fes)
54-
55-
# end
56-
# end
57-
5819
# TODO
5920
# 1. need to loop over bcs and vars in dof
6021
# to make organize different bcs based on fspace type
6122
# 2. group all dbcs, nbcs of similar fspaces into a single struct?
62-
# 3. figure out how to handle function pointers on the GPU
23+
# 3. figure out how to handle function pointers on the GPU - done
6324
# 4. add different fspace types
64-
# 5. convert vectors of dbcs/nbcs into namedtuples
25+
# 5. convert vectors of dbcs/nbcs into namedtuples - done
6526
function Parameters(
6627
assembler, physics,
6728
dirichlet_bcs,
@@ -76,6 +37,16 @@ function Parameters(
7637
state_old = nothing
7738
state_new = nothing
7839

40+
# for mixed spaces we'll need to do this more carefully
41+
if isa(physics, AbstractPhysics)
42+
syms = keys(values(assembler.dof.H1_vars)[1].fspace.elem_conns)
43+
physics = map(x -> physics, syms)
44+
physics = NamedTuple{tuple(syms...)}(tuple(physics...))
45+
else
46+
@assert isa(physics, NamedTuple)
47+
# TODO re-arrange physics tuple to match fspaces when appropriate
48+
end
49+
7950
if dirichlet_bcs !== nothing
8051
syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dirichlet_bcs))
8152
# dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))

‎src/TimeSteppers.jl‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
abstract type AbstractTimeStepper{T} end
2+
current_time(t::AbstractTimeStepper) = sum(t.time_current)
3+
time_step(t::AbstractTimeStepper) = sum(t.Δt)
4+
5+
struct TimeStepper{T} <: AbstractTimeStepper{T}
6+
time_start::T
7+
time_end::T
8+
time_current::T
9+
Δt::T
10+
end
11+
12+
function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
13+
time_start = zeros(1)
14+
time_end = zeros(1)
15+
time_current = zeros(1)
16+
Δt = zeros(1)
17+
Δt = zeros(1)
18+
fill!(time_start, time_start_in)
19+
fill!(time_end, time_end_in)
20+
fill!(time_current, time_start_in)
21+
fill!(Δt, (time_end_in - time_start_in) / n_steps)
22+
return TimeStepper(time_start, time_end, time_current, Δt)
23+
# return TimeStepper(time_start_in, time_end_in, )
24+
end

0 commit comments

Comments
 (0)