Inclined Support for Single Quasi-Static Problems

examples/single/static-solid/cube_compression_inclined_support_30_degrees/cube.jou

@@ -0,0 +1,51 @@
+${side = 1.0}
+${length = 1.0}
+${separation = 1.0}
+${nel_fine = 2}
+${nel_coarse = 1}
+${h_fine = length / nel_fine}
+${h_coarse = length / nel_coarse}
+#{cos_th = 0.8660254037844387}
+#{sin_th = .5}
+${area = side * side}
+${offset_x = (length + separation)*cos_th / 2.0}
+${offset_y = (length + separation)*sin_th / 2.0}
+
+create brick x {length} y {side} z {side}
+move volume 1 x -0.5 include_merged
+move volume 1 y -0.5 include_merged
+#volume 1 scheme tetmesh
+volume 1 size {h_fine}
+mesh volume 1
+block 1 volume 1
+#block 1 element type tetra4
+block 1 name "fine"
+nodeset 1 surface 4
+nodeset 1 name "nsx-"
+nodeset 2 surface 6
+nodeset 2 name "nsx+"
+nodeset 3 surface 3
+nodeset 3 name "nsy-"
+nodeset 4 surface 5
+nodeset 4 name "nsy+"
+nodeset 5 surface 2
+nodeset 5 name "nsz-"
+nodeset 6 surface 1
+nodeset 6 name "nsz+"
+nodeset 7 volume all
+nodeset 7 name "nsall"
+sideset 1 surface 4
+sideset 1 name "ssx-"
+sideset 2 surface 6
+sideset 2 name "ssx+"
+sideset 3 surface 3
+sideset 3 name "ssy-"
+sideset 4 surface 5
+sideset 4 name "ssy+"
+sideset 5 surface 2
+sideset 5 name "ssz-"
+sideset 6 surface 1
+sideset 6 name "ssz+"
+set large exodus file off
+rotate volume 1 about z angle 30
+export mesh "cube-1.g" overwrite

examples/single/static-solid/cube_compression_inclined_support_30_degrees/cube.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:
+      fine: 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
+  final time: 2
+  time step: 0.5
+boundary conditions:
+  Dirichlet:
+    - node set: nsx-
+      component: x
+      function: "0.866025403784439 * 0.1*t"
+    - node set: nsx-
+      component: y
+      function: "0.5 * 0.1*t"
+    - node set: nsx-
+      component: z
+      function: "0"
+  Inclined Dirichlet:
+    - node set: nsx+
+      normal vector: [ 0.866025403784439, 0.5, 0 ]
+      function: "-0.1*t"
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 50
+  relative tolerance: 1.0e-12
+  absolute tolerance: 1.0e-08

src/ics_bcs.jl

@@ -22,6 +22,42 @@ function SMDirichletBC(input_mesh::ExodusDatabase, bc_params::Dict{Any,Any})
     )
 end
 
+function SMDirichletInclined(input_mesh::ExodusDatabase, bc_params::Dict{Any,Any})
+    node_set_name = bc_params["node set"]
+    expression = bc_params["function"]
+    node_set_id = node_set_id_from_name(node_set_name, input_mesh)
+    node_set_node_indices = Exodus.read_node_set_nodes(input_mesh, node_set_id)
+    # expression is an arbitrary function of t, x, y, z in the input file
+    disp_num = eval(Meta.parse(expression))
+    velo_num = expand_derivatives(D(disp_num))
+    acce_num = expand_derivatives(D(velo_num))
+    # For inclined support, the function is applied along the x direction
+    offset = component_offset_from_string("x")
+
+    # The local basis is determined from a normal vector
+    axis = bc_params["normal vector"]
+    axis = axis/norm(axis)
+    e1 = [1.0, 0.0, 0.0]
+    w = cross(e1, axis)
+    s = norm(w)
+    θ = asin(s)
+    m = w/s
+    rv = θ * m
+    # Rotation is converted via the psuedo vector to rotation matrix
+    rotation_matrix = MiniTensor.rt_from_rv(rv)'
+
+    SMDirichletInclined(
+        node_set_name,
+        node_set_id,
+        node_set_node_indices,
+        disp_num,
+        velo_num,
+        acce_num,
+        rotation_matrix,
+        offset
+    )
+end
+
 function SMNeumannBC(input_mesh::ExodusDatabase, bc_params::Dict{Any,Any})
     side_set_name = bc_params["side set"]
     expression = bc_params["function"]
@@ -183,6 +219,43 @@ function apply_bc(model::SolidMechanics, bc::SMNeumannBC)
     end
 end
 
+function apply_bc(model::SolidMechanics, bc::SMDirichletInclined)
+    for node_index ∈ bc.node_set_node_indices
+        values = Dict(
+            t => model.time,
+            x => model.reference[1, node_index],
+            y => model.reference[2, node_index],
+            z => model.reference[3, node_index],
+        )
+        disp_sym = substitute(bc.disp_num, values)
+        velo_sym = substitute(bc.velo_num, values)
+        acce_sym = substitute(bc.acce_num, values)
+
+        # For inclined support, these values are in the local direction
+        disp_val_loc = extract_value(disp_sym)
+        velo_val_loc = extract_value(velo_sym)
+        acce_val_loc = extract_value(acce_sym)
+
+        # Rotate the local displacements to the global frame
+        disp_vector_local = zeros(3)
+        disp_vector_local[bc.offset] = disp_val_loc 
+        velo_vector_local = zeros(3)
+        velo_vector_local[bc.offset] = velo_val_loc 
+        accel_vector_local = zeros(3)
+        accel_vector_local[bc.offset] = acce_val_loc     
+        disp_vector_glob = bc.rotation_matrix' * disp_vector_local
+        velo_vector_glob = bc.rotation_matrix' * velo_vector_local
+        accel_vector_glob = bc.rotation_matrix' * accel_vector_local
+
+        dof_index = 3 * (node_index - 1) + bc.offset
+        model.current[:, node_index] =
+            model.reference[:, node_index] + disp_vector_glob
+        model.velocity[:, node_index] = velo_vector_glob
+        model.acceleration[:, node_index] = accel_vector_glob
+        model.free_dofs[dof_index] = false
+    end
+end
+
 function find_point_in_mesh(
     point::Vector{Float64},
     model::SolidMechanics,
@@ -541,6 +614,7 @@ function create_bcs(params::Dict{Any,Any})
     end
     input_mesh = params["input_mesh"]
     bc_params = params["boundary conditions"]
+    inclined_support_nodes = Vector{Int64}()
     for (bc_type, bc_type_params) ∈ bc_params
         for bc_setting_params ∈ bc_type_params
             if bc_type == "Dirichlet"
@@ -558,6 +632,10 @@ function create_bcs(params::Dict{Any,Any})
                 boundary_condition =
                     SMContactSchwarzBC(coupled_subsim, input_mesh, bc_setting_params)
                 push!(boundary_conditions, boundary_condition)
+            elseif bc_type == "Inclined Dirichlet"
+                boundary_condition = SMDirichletInclined(input_mesh, bc_setting_params)
+                append!(inclined_support_nodes, boundary_condition.node_set_node_indices)
+                push!(boundary_conditions, boundary_condition)
             elseif bc_type == "Schwarz overlap" || bc_type == "Schwarz nonoverlap"
                 sim = params["global_simulation"]
                 subsim_name = params["name"]
@@ -574,6 +652,11 @@ function create_bcs(params::Dict{Any,Any})
             end
         end
     end
+    # BRP: do not support applying multiple inclined support BCs to a single node
+    duplicate_inclined_support_conditions = length(unique(inclined_support_nodes)) < length(inclined_support_nodes)
+    if duplicate_inclined_support_conditions
+        throw(error("Cannot apply multiple inclined BCs to a single node."))
+    end
     return boundary_conditions
 end
 

src/ics_bcs_def.jl

@@ -20,6 +20,17 @@ mutable struct SMDirichletBC <: RegularBoundaryCondition
     acce_num::Num
 end
 
+mutable struct SMDirichletInclined <: RegularBoundaryCondition
+    node_set_name::String
+    node_set_id::Int64
+    node_set_node_indices::Vector{Int64}
+    disp_num::Num
+    velo_num::Num
+    acce_num::Num
+    rotation_matrix::Matrix{Float64}
+    offset::Int64
+end
+
 mutable struct SMNeumannBC <: RegularBoundaryCondition
     side_set_name::String
     offset::Int64

src/minitensor.jl

@@ -1,4 +1,5 @@
 module MiniTensor
+using LinearAlgebra: norm
 
 #
 # Lie groups and Lie algebra utilities, mostly algebra of rotations

src/model.jl

@@ -76,6 +76,11 @@ function SolidMechanics(params::Dict{Any,Any})
     else
         smooth_reference = ""
     end
+
+    # BRP: define a global transform for inclined support
+    global_transform_t = Diagonal(ones(3 * num_nodes))
+    global_transform = sparse(global_transform_t)
+
     SolidMechanics(
         input_mesh,
         materials,
@@ -93,6 +98,7 @@ function SolidMechanics(params::Dict{Any,Any})
         failed,
         mesh_smoothing,
         smooth_reference,
+        global_transform
     )
 end
 
@@ -308,6 +314,19 @@ function evaluate(integrator::QuasiStatic, model::SolidMechanics)
     stiffness = Vector{Float64}()
     blocks = Exodus.read_sets(input_mesh, Block)
     num_blks = length(blocks)
+
+    # BRP: bookkeeping parameters for inclined support BCs
+    inclined_support_node_indices = Vector{Int64}()
+    inclined_support_bc_indices = Vector{Int64}()
+    for (bc_idx, bc) in enumerate(model.boundary_conditions)
+        if isa(bc, SMDirichletInclined)
+            append!(inclined_support_node_indices, bc.node_set_node_indices)
+            for _ in bc.node_set_node_indices
+                push!(inclined_support_bc_indices, bc_idx)
+            end
+        end
+    end
+
     for blk_index ∈ 1:num_blks
         material = materials[blk_index]
         ρ = material.ρ
@@ -367,6 +386,16 @@ function evaluate(integrator::QuasiStatic, model::SolidMechanics)
             assemble(rows, cols, stiffness, element_stiffness, elem_dofs)
         end
     end
+
+    # BRP Inclined Support Operations
+    T_local = Matrix(Diagonal(ones(num_dof)))
+    for (corresponding_bc_idx, inc_support_node_idx) in zip(inclined_support_bc_indices, inclined_support_node_indices)
+        T_nodal = model.boundary_conditions[corresponding_bc_idx].rotation_matrix
+        base = 3*(inc_support_node_idx-1) # Block index in global stiffness
+        T_local[base+1:base+3, base+1:base+3] *= T_nodal
+    end
+    model.global_transform = sparse(T_local)
+
     stiffness_matrix = sparse(rows, cols, stiffness)
     model.internal_force = internal_force
     return energy, internal_force, body_force, stiffness_matrix

src/model_def.jl

@@ -1,4 +1,5 @@
 abstract type Model end
+using SparseArrays
 
 mutable struct SolidMechanics <: Model
     mesh::ExodusDatabase
@@ -17,6 +18,7 @@ mutable struct SolidMechanics <: Model
     failed::Bool
     mesh_smoothing::Bool
     smooth_reference::String
+    global_transform::SparseArrays.Matrix{Float64}
 end
 
 # TODO: Add potential energy as in the above

src/solver.jl

@@ -182,8 +182,11 @@ function copy_solution_source_targets(
     solver::Any,
     model::SolidMechanics,
 )
-    displacement = integrator.displacement
-    solver.solution = displacement
+    displacement_local = integrator.displacement
+    solver.solution = displacement_local
+    # BRP: apply inclined support inverse transform
+    displacement = model.global_transform' * displacement_local
+
     _, num_nodes = size(model.reference)
     for node ∈ 1:num_nodes
         nodal_displacement = displacement[3*node-2:3*node]
@@ -196,8 +199,9 @@ function copy_solution_source_targets(
     model::SolidMechanics,
     integrator::QuasiStatic,
 )
-    displacement = solver.solution
-    integrator.displacement = displacement
+    displacement_local = solver.solution
+    integrator.displacement = displacement_local
+    displacement = model.global_transform' * displacement_local
     _, num_nodes = size(model.reference)
     for node ∈ 1:num_nodes
         nodal_displacement = displacement[3*node-2:3*node]
@@ -215,6 +219,8 @@ function copy_solution_source_targets(
         nodal_displacement = model.current[:, node] - model.reference[:, node]
         integrator.displacement[3*node-2:3*node] = nodal_displacement
     end
+    # Convert integrator displacement from global to local
+    integrator.displacement = model.global_transform * integrator.displacement
     solver.solution = integrator.displacement
 end
 
@@ -341,8 +347,8 @@ function evaluate(integrator::QuasiStatic, solver::HessianMinimizer, model::Soli
     integrator.stored_energy = stored_energy
     solver.value = stored_energy
     external_force = body_force + model.boundary_force
-    solver.gradient = internal_force - external_force
-    solver.hessian = stiffness_matrix
+    solver.gradient = model.global_transform * (internal_force - external_force)
+    solver.hessian = model.global_transform * stiffness_matrix
 end
 
 function evaluate(integrator::QuasiStatic, solver::SteepestDescent, model::SolidMechanics)
@@ -577,4 +583,5 @@ function solve(integrator::TimeIntegrator, solver::Solver, model::Model)
             break
         end
     end
+    solver.gradient = model.global_transform' * solver.gradient
 end

test/runtests.jl

@@ -19,3 +19,4 @@ include("tet10-static-solid-cube.jl")
 include("schwarz-overlap-static-cuboid-hex8.jl")
 include("schwarz-nonoverlap-static-cuboid-hex8.jl")
 include("schwarz-contact-static-cubes-hex8.jl")
+include("single-static-solid-cube-inclined-support.jl")