Physics/block by block

examples/mechanics/Cube.jl

@@ -0,0 +1,88 @@
+using Exodus
+using FiniteElementContainers
+using Tensors
+
+fixed(_, _) = 0.
+displace(_, t) = 1.0 * t
+
+struct Mechanics{Form} <: AbstractPhysics{2, 0, 0}
+  formulation::Form
+end
+
+function strain_energy(∇u)
+  E = 1.e9
+  ν = 0.25
+  K = E / (3. * (1. - 2. * ν))
+  G = E / (2. * (1. + ν))
+
+  ε = symmetric(∇u)
+  ψ = 0.5 * K * tr(ε)^2 + G * dcontract(dev(ε), dev(ε))
+end
+
+function FiniteElementContainers.residual(physics::Mechanics, cell, u_el, args...)
+  (; X_q, N, ∇N_X, JxW) = cell
+
+  # kinematics
+  ∇u_q = u_el * ∇N_X
+  ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
+  # constitutive
+  P_q = gradient(strain_energy, ∇u_q)
+  # turn into voigt notation
+  P_q = extract_stress(physics.formulation, P_q)
+  G_q = discrete_gradient(physics.formulation, ∇N_X)
+  f_q = G_q * P_q
+  return JxW * f_q[:]
+  # return f_q
+end
+
+function FiniteElementContainers.stiffness(physics::Mechanics, cell, u_el, args...)
+  (; X_q, N, ∇N_X, JxW) = cell
+
+  # kinematics
+  ∇u_q = u_el * ∇N_X
+  ∇u_q = modify_field_gradients(physics.formulation, ∇u_q)
+  # constitutive
+  K_q = hessian(z -> strain_energy(z), ∇u_q)
+  # turn into voigt notation
+  K_q = extract_stiffness(physics.formulation, K_q)
+  G_q = discrete_gradient(physics.formulation, ∇N_X)
+  return JxW * G_q * K_q * G_q'
+end
+
+function mechanics_test()
+  mesh = UnstructuredMesh("cube.g")
+  V = FunctionSpace(mesh, H1Field, Lagrange) 
+  physics = (;
+    cube=Mechanics(ThreeDimensional())
+  )
+
+  u = VectorFunction(V, :displ)
+  asm = SparseMatrixAssembler(H1Field, u)
+
+  dbcs = DirichletBC[
+    DirichletBC("displ_x", "ssx-", fixed),
+    DirichletBC("displ_y", "ssy-", fixed),
+    DirichletBC("displ_z", "ssz-", fixed),
+    DirichletBC("displ_z", "ssz+", displace),
+  ]
+
+  # pre-setup some scratch arrays
+  times = TimeStepper(0., 1., 10)
+  p = create_parameters(asm, physics; dirichlet_bcs=dbcs, times=times)
+
+  solver = NewtonSolver(IterativeLinearSolver(asm, :CgSolver))
+  integrator = QuasiStaticIntegrator(solver)
+  pp = PostProcessor(mesh, "cube.e", u)
+
+  write_times(pp, 1, p.times.time_current[1])
+  write_field(pp, 1, p.h1_field)
+
+  for n in 1:10
+    evolve!(integrator, p)
+    write_times(pp, n, p.times.time_current[1])
+    write_field(pp, n, p.h1_field)
+  end
+  close(pp)
+end
+
+mechanics_test()

examples/mechanics/cube-line-search.yaml

@@ -0,0 +1,45 @@
+type: single
+input mesh file: cube.g
+output mesh file: cube.e
+Exodus output interval: 1
+CSV output interval: 0
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      cube: elastic
+    elastic:
+      model: linear elastic
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: quasi static
+  initial time: 0.0
+  final time: 1.0
+  time step: 0.1
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  use line search: true
+  line search backtrack factor: 0.5
+  line search decrease factor: 1.0e-04
+  line search maximum iterations: 16
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-12
+  absolute tolerance: 1.0e-08

examples/mechanics/cube-steepest-descent.yaml

@@ -0,0 +1,42 @@
+type: single
+input mesh file: cube.g
+output mesh file: cube.e
+Exodus output interval: 1
+CSV output interval: 0
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      cube: elastic
+    elastic:
+      model: linear elastic
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: quasi static
+  initial time: 0.0
+  final time: 1.0
+  time step: 0.1
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+solver:
+  type: steepest descent
+  step: steepest descent
+  step length: 1.0
+  minimum iterations: 1
+  maximum iterations: 1024
+  relative tolerance: 1.0e-15
+  absolute tolerance: 1.0e-08

examples/mechanics/cube.yaml

@@ -0,0 +1,41 @@
+type: single
+input mesh file: cube.g
+output mesh file: cube.e
+Exodus output interval: 1
+CSV output interval: 0
+model:
+  type: solid mechanics
+  material:
+    blocks:
+      cube: elastic
+    elastic:
+      model: linear elastic
+      elastic modulus: 1.0e+09
+      Poisson's ratio: 0.25
+      density: 1000.0
+time integrator:
+  type: quasi static
+  initial time: 0.0
+  final time: 1.0
+  time step: 0.1
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.0"
+    - node set: nsy-
+      component: y
+      function: "0.0"
+    - node set: nsz-
+      component: z
+      function: "0.0"
+    - node set: nsz+
+      component: z
+      function: "1.0 * t"
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-12
+  absolute tolerance: 1.0e-08

src/DofManagers.jl

@@ -219,6 +219,8 @@ Base.length(dof::DofManager) = length(dof.H1_bc_dofs) + length(dof.H1_unknown_do
 
 KA.get_backend(dof::DofManager) = KA.get_backend(dof.H1_unknown_dofs)
 
+# is create_bcs even really needed anymore?
+
 """
 $(TYPEDSIGNATURES)
 """

src/FiniteElementContainers.jl

@@ -143,6 +143,7 @@ include("physics/Physics.jl")
 include("PostProcessors.jl")
 
 #
+include("TimeSteppers.jl")
 include("Parameters.jl")
 include("solvers/Solvers.jl")
 include("integrators/Integrators.jl")

src/Parameters.jl

@@ -1,29 +1,3 @@
-# move to a seperate file
-abstract type AbstractTimeStepper{T} end
-current_time(t::AbstractTimeStepper) = sum(t.time_current)
-time_step(t::AbstractTimeStepper) = sum(t.Δt)
-
-struct TimeStepper{T} <: AbstractTimeStepper{T}
-  time_start::T
-  time_end::T
-  time_current::T
-  Δt::T
-end
-
-function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
-  time_start = zeros(1)
-  time_end = zeros(1)
-  time_current = zeros(1)
-  Δt = zeros(1)
-  Δt = zeros(1)
-  fill!(time_start, time_start_in)
-  fill!(time_end, time_end_in)
-  fill!(time_current, time_start_in)
-  fill!(Δt, (time_end_in - time_start_in) / n_steps)
-  return TimeStepper(time_start, time_end, time_current, Δt)
-  # return TimeStepper(time_start_in, time_end_in, )
-end
-
 abstract type AbstractParameters end
 
 struct Parameters{D, N, T, Phys, Props, S, V, H1} <: AbstractParameters
@@ -42,26 +16,13 @@
   h1_field::H1
 end
 
-# TODO only works for H1Fields currently most likely
-# function Parameters(dof::DofManager, physics::AbstractPhysics)
-#   n_props = num_properties(physics)
-#   n_state = num_states(physics)
-
-#   # TODO
-#   fspace = dof.H1_vars[1].fspace
-
-#   for (name, ref_fe) in pairs(fspace.ref_fes)
-
-#   end
-# end
-
 # TODO 
 # 1. need to loop over bcs and vars in dof
 #    to make organize different bcs based on fspace type
 # 2. group all dbcs, nbcs of similar fspaces into a single struct?
-# 3. figure out how to handle function pointers on the GPU
+# 3. figure out how to handle function pointers on the GPU - done
 # 4. add different fspace types
-# 5. convert vectors of dbcs/nbcs into namedtuples
+# 5. convert vectors of dbcs/nbcs into namedtuples - done
 function Parameters(
   assembler, physics,
   dirichlet_bcs, 
@@ -76,6 +37,16 @@
   state_old = nothing
   state_new = nothing
 
+  # for mixed spaces we'll need to do this more carefully
+  if isa(physics, AbstractPhysics)
+    syms = keys(values(assembler.dof.H1_vars)[1].fspace.elem_conns)
+    physics = map(x -> physics, syms)
+    physics = NamedTuple{tuple(syms...)}(tuple(physics...))
+  else
+    @assert isa(physics, NamedTuple)
+    # TODO re-arrange physics tuple to match fspaces when appropriate
+  end
+
   if dirichlet_bcs !== nothing
     syms = map(x -> Symbol("dirichlet_bc_$x"), 1:length(dirichlet_bcs))
     # dbcs = NamedTuple{tuple(syms...)}(tuple(dbcs...))

src/TimeSteppers.jl

@@ -0,0 +1,24 @@
+abstract type AbstractTimeStepper{T} end
+current_time(t::AbstractTimeStepper) = sum(t.time_current)
+time_step(t::AbstractTimeStepper) = sum(t.Δt)
+
+struct TimeStepper{T} <: AbstractTimeStepper{T}
+  time_start::T
+  time_end::T
+  time_current::T
+  Δt::T
+end
+
+function TimeStepper(time_start_in::T, time_end_in::T, n_steps::Int) where T <: Number
+  time_start = zeros(1)
+  time_end = zeros(1)
+  time_current = zeros(1)
+  Δt = zeros(1)
+  Δt = zeros(1)
+  fill!(time_start, time_start_in)
+  fill!(time_end, time_end_in)
+  fill!(time_current, time_start_in)
+  fill!(Δt, (time_end_in - time_start_in) / n_steps)
+  return TimeStepper(time_start, time_end, time_current, Δt)
+  # return TimeStepper(time_start_in, time_end_in, )
+end