Merge pull request #128 from sandialabs/pressure_bc

Adding Pressure NBC to Norma

Author: lxmota

examples/single/static-solid/pressure-bc/cube.g

@@ -0,0 +1,34 @@
+${side = 1.0}
+${h_size = 0.5}
+
+create brick x {side} y {side} z {side}
+volume all size {h_size}
+mesh volume all
+block 1 volume 1
+block 1 name "cube"
+nodeset 1 surface with x_coord < 0
+nodeset 1 name "nsx-"
+nodeset 2 surface with x_coord > 0
+nodeset 2 name "nsx+"
+nodeset 3 surface with y_coord < 0
+nodeset 3 name "nsy-"
+nodeset 4 surface with y_coord > 0
+nodeset 4 name "nsy+"
+nodeset 5 surface with z_coord < 0
+nodeset 5 name "nsz-"
+nodeset 6 surface with z_coord > 0
+nodeset 6 name "nsz+"
+sideset 1 surface with x_coord < 0
+sideset 1 name "ssx-"
+sideset 2 surface with x_coord > 0
+sideset 2 name "ssx+"
+sideset 3 surface with y_coord < 0
+sideset 3 name "ssy-"
+sideset 4 surface with y_coord < 0
+sideset 4 name "ssy+"
+sideset 5 surface with z_coord < 0
+sideset 5 name "ssz-"
+sideset 6 surface with z_coord > 0
+sideset 6 name "ssz+"
+set large exodus file off
+export mesh "cube.g" overwrite

examples/single/static-solid/pressure-bc/cube.jou

@@ -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"
+  Neumann pressure:
+    - side set: ssz+
+      function: "1.0e+09 * t"
+solver:
+  type: Hessian minimizer
+  step: full Newton
+  minimum iterations: 1
+  maximum iterations: 16
+  relative tolerance: 1.0e-14
+  absolute tolerance: 1.0e-10

examples/single/static-solid/pressure-bc/cube.yaml

@@ -77,6 +77,24 @@ function SolidMechanicsNeumannBoundaryCondition(input_mesh::ExodusDatabase, bc_p
     )
 end
 
+function SolidMechanicsNeumannPressureBoundaryCondition(input_mesh::ExodusDatabase, bc_params::Parameters)
+    side_set_name = bc_params["side set"]
+    expression = bc_params["function"]
+    side_set_id = side_set_id_from_name(side_set_name, input_mesh)
+    num_nodes_per_side, side_set_node_indices = Exodus.read_side_set_node_list(input_mesh, side_set_id)
+    side_set_node_indices = Int64.(side_set_node_indices)
+
+    # Build symbolic expressions
+    pressure_num = eval(Meta.parse(expression))
+
+    # Compile them into functions
+    pressure_fun = eval(build_function(pressure_num, [t, x, y, z]; expression=Val(false)))
+
+    return SolidMechanicsNeumannPressureBoundaryCondition(
+        side_set_name, side_set_id, num_nodes_per_side, side_set_node_indices, pressure_fun
+    )
+end 
+
 function SolidMechanicsOverlapSchwarzBoundaryCondition(
     coupled_block_name::String,
     tol::Float64,
@@ -313,6 +331,22 @@ function apply_bc(model::SolidMechanics, bc::SolidMechanicsNeumannBoundaryCondit
     end
 end
 
+function apply_bc(model::SolidMechanics, bc::SolidMechanicsNeumannPressureBoundaryCondition)
+    ss_node_index = 1
+    for side in bc.num_nodes_per_side
+        side_nodes = bc.side_set_node_indices[ss_node_index:(ss_node_index + side - 1)]
+        side_coordinates = model.reference[:, side_nodes]
+        nodal_force_component = get_side_set_nodal_pressure(side_coordinates, bc.pressure_fun, model.time)
+        ss_node_index += side
+        side_node_index = 1
+        for node_index in side_nodes
+            dof_indices = [3 * node_index - 2, 3 * node_index - 1,  3 * node_index]
+            model.boundary_force[dof_indices] += nodal_force_component[:, side_node_index]
+            side_node_index += 1
+        end
+    end
+end
+
 function compute_rotation_matrix(axis::SVector{3,Float64})::SMatrix{3,3,Float64}
     e1 = @SVector [1.0, 0.0, 0.0]
     angle_btwn = acos(dot(axis, e1))
@@ -846,6 +880,9 @@ function create_bcs(params::Parameters)
             elseif bc_type == "Neumann"
                 boundary_condition = SolidMechanicsNeumannBoundaryCondition(input_mesh, bc_setting_params)
                 push!(boundary_conditions, boundary_condition)
+            elseif bc_type == "Neumann pressure"
+                boundary_condition = SolidMechanicsNeumannPressureBoundaryCondition(input_mesh, bc_setting_params)
+                push!(boundary_conditions, boundary_condition)
             elseif bc_type == "Schwarz contact"
                 sim = params["parent_simulation"]
                 sim.controller.schwarz_contact = true

src/ics_bcs.jl

@@ -42,6 +42,17 @@ mutable struct SolidMechanicsNeumannBoundaryCondition <: SolidMechanicsRegularBo
     traction_fun::Function
 end
 
+#IKT 6/9/2025 TODO: check with Alejandro if want to have separate NeumannPressure struct
+#or integrate it into the Neumann struct.  The latter is possible and avoids some
+#code duplication.
+mutable struct SolidMechanicsNeumannPressureBoundaryCondition <: SolidMechanicsRegularBoundaryCondition
+    name::String
+    side_set_id::Int64
+    num_nodes_per_side::Vector{Int64}
+    side_set_node_indices::Vector{Int64}
+    pressure_fun::Function
+end
+
 mutable struct SolidMechanicsContactSchwarzBoundaryCondition <: SolidMechanicsSchwarzBoundaryCondition
     name::String
     side_set_id::Int64

src/ics_bcs_types.jl

@@ -358,6 +358,29 @@ function get_side_set_nodal_forces(nodal_coord::Matrix{Float64}, traction_fun::F
     return nodal_force_component
 end
 
+function get_side_set_nodal_pressure(nodal_coord::Matrix{Float64}, pressure_fun::Function, time::Float64)
+    _, num_side_nodes = size(nodal_coord)
+    element_type = get_element_type(2, num_side_nodes)
+    num_int_points = default_num_int_pts(element_type)
+    N, dNdξ, w, _ = isoparametric(element_type, num_int_points)
+    nodal_force_component = zeros(3, num_side_nodes)
+    for point in 1:num_int_points
+        Nₚ = N[:, point]
+        dNdξₚ = dNdξ[:, :, point]
+        dXdξ = dNdξₚ * nodal_coord'
+        perp_vector = cross(dXdξ[1, :], dXdξ[2, :])
+        normal = LinearAlgebra.normalize(perp_vector)  
+        j = norm(perp_vector)
+        wₚ = w[point]
+        point_coord = nodal_coord * Nₚ
+        txzy = (time, point_coord[1], point_coord[2], point_coord[3])
+        pressure_val = pressure_fun(txzy...)
+        nodal_force_component_vector = LinearAlgebra.kron(pressure_val * Nₚ * j * wₚ, normal)
+        nodal_force_component += reshape(nodal_force_component_vector, (3, num_side_nodes))
+    end
+    return nodal_force_component
+end
+
 function map_to_parametric(element_type::ElementType, nodes::Matrix{Float64}, point::Vector{Float64})
     tol = 1.0e-08
     max_iters = 1024

src/interpolation.jl

@@ -55,7 +55,8 @@ const indexed_test_files = [
     (39, "schwarz-ahead-nonoverlap-dynamic-torsion.jl"),
     (40, "schwarz-ahead-nonoverlap-dynamic-plate.jl"),
     (41, "schwarz-ahead-nonoverlap-dynamic-bracket.jl"),
-    (42, "utils.jl"), # Must go last due to FPE traps
+    (42, "single-static-solid-pressure-bc.jl"), 
+    (43, "utils.jl"), # Must go last due to FPE traps
 ]
 
 # Extract test file names

test/runtests.jl

@@ -0,0 +1,26 @@
+# Norma: Copyright 2025 National Technology & Engineering Solutions of
+# Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
+# the U.S. Government retains certain rights in this software. This software
+# is released under the BSD license detailed in the file license.txt in the
+# top-level Norma.jl directory.
+@testset "Single Static Solid Pressure Neumann Bc" begin
+    cp("../examples/single/static-solid/pressure-bc/cube.yaml", "cube.yaml"; force=true)
+    cp("../examples/single/static-solid/pressure-bc/cube.g", "cube.g"; force=true)
+    simulation = Norma.run("cube.yaml")
+    integrator = simulation.integrator
+    model = simulation.model
+    rm("cube.yaml")
+    rm("cube.g")
+    rm("cube.e")
+    avg_disp = average_components(integrator.displacement)
+    avg_stress = average_components(model.stress)
+    @test avg_disp[1] ≈ -0.125 rtol = 1.0e-06
+    @test avg_disp[2] ≈ -0.125 rtol = 1.0e-06
+    @test avg_disp[3] ≈ 0.500 rtol = 1.0e-06
+    @test avg_stress[1] ≈ 0.0 atol = 1.0e-06
+    @test avg_stress[2] ≈ 0.0 atol = 1.0e-06
+    @test avg_stress[3] ≈ 1.0e+09 rtol = 1.0e-06
+    @test avg_stress[4] ≈ 0.0 atol = 1.0e-06
+    @test avg_stress[5] ≈ 0.0 atol = 1.0e-06
+    @test avg_stress[6] ≈ 0.0 atol = 1.0e-06
+end