Skip to content

add python mesh functions using gmsh. This will be used to create in-built meshing capabilities in PeriDEM for specific geometries. #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
prashjha merged 1 commit into from
Jun 3, 2025
Merged

add python mesh functions using gmsh. This will be used to create in-built meshing capabilities in PeriDEM for specific geometries. #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
273 changes: 273 additions & 0 deletions tools/python_utils/circle_wall_input_using_symmetric_mesh.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,273 @@
import os
import sys
import numpy as np

#sys.path.append('<Path to PeriDEM source>/tools/python_utils')
from gmsh_particles import *
from util import *

def create_input_file(inp_dir, pp_tag):
"""Generates input file for two-particle test"""

sim_inp_dir = str(inp_dir)

R1 = 0.001
mesh_size = R1 / 8.
horizon = 2.2 * mesh_size
particle_dist = 0.9*mesh_size #horizon # surface to surface distance

## particle 1 rectangle
Lx, Ly = 4*R1, 3*mesh_size
rect_center = [R1, 0., 0.]

## particle 2 circle
cir_center = [R1, rect_center[1] + 0.5*Ly + particle_dist + R1, 0.]

## assign free fall velocity to second particle
high_impact = False
free_fall_vel = [0., -0.1, 0.]
if high_impact:
free_fall_vel = [0., -4., 0.] # high impact velocity

## time
final_time = 0.01
num_steps = 50000
if high_impact:
final_time = 0.0002
num_steps = 10000
# final_time = 0.00002
# num_steps = 2
num_outputs = 100
dt_out_n = num_steps / num_outputs
perform_out = True

## material
poisson1 = 0.25
rho1 = 1200.
K1 = 2.16e+7
E1 = get_E(K1, poisson1)
G1 = get_G(E1, poisson1)
Gc1 = 50.

poisson2 = 0.25
rho2 = 1200.
K2 = 2.16e+7
E2 = get_E(K2, poisson2)
G2 = get_G(E2, poisson2)
Gc2 = 50.

## contact
# R_contact = 0.95 * mesh_size
# R_contact = 1.74e-04
R_contact_factor = 0.95

# Kn_V_max = 7.385158e+05
# Kn = np.power(Kn_V_max, 2)
# compute from bulk modulus

# from bulk modulus
Kn_11 = 18. * get_eff_k(K1, K1) / (np.pi * np.power(horizon, 5))
Kn_22 = 18. * get_eff_k(K2, K2) / (np.pi * np.power(horizon, 5))
Kn_12 = 18. * get_eff_k(K1, K2) / (np.pi * np.power(horizon, 5))

beta_n_eps = 0.9
friction_coeff = 0.5
damping_active = False
friction_active = False
beta_n_factor = 10.
Kn_factor = 0.1

## gravity
gravity_active = True
gravity = [0., -10., 0.]

## neighbor search details
neigh_search_factor = 2.
neigh_search_interval = 1
neigh_search_criteria = "simple_all"

### ---------------------------------------------------------------- ###
# generate mesh and particle location data
### ---------------------------------------------------------------- ###
plocf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
plocf.write("i, x, y, z, r, o\n")
plocf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (0, rect_center[0], rect_center[1], rect_center[2], Lx, 0.))
plocf.write("%d, %Lf, %Lf, %Lf, %Lf, %Lf\n" % (1, cir_center[0], cir_center[1], cir_center[2], R1, 0.))
plocf.close()

zones_mesh_fnames = ["mesh_rect_1", "mesh_cir_2"]

# generate mesh for particle 1
rectangle_mesh_symmetric(xc = [0., 0., 0.], Lx = Lx, Ly = Ly, h = mesh_size, filename = zones_mesh_fnames[0] + "_" + str(pp_tag), vtk_out = True, symmetric_mesh = True)

# generate mesh for particle 2
circle_mesh_symmetric(xc = [0., 0., 0.], r = R1, h = mesh_size, filename = zones_mesh_fnames[1] + "_" + str(pp_tag), vtk_out = True, symmetric_mesh = True)

os.system("mkdir -p ../out")

### ---------------------------------------------------------------- ###
# generate YAML file
### ---------------------------------------------------------------- ###
# print('\nGenerating imput file\n')
inpf = open(sim_inp_dir + 'input_' + str(pp_tag) + '.yaml','w')
inpf.write("Model:\n")
inpf.write(" Dimension: 2\n")
inpf.write(" Discretization_Type:\n")
inpf.write(" Spatial: finite_difference\n")
inpf.write(" Time: central_difference\n")
inpf.write(" Final_Time: %4.6e\n" % (final_time))
inpf.write(" Time_Steps: %d\n" % (num_steps))

#
# container info
#
inpf.write("Container:\n")
inpf.write(" Geometry:\n")
inpf.write(" Type: rectangle\n")
contain_params = [rect_center[0] - 0.5*Lx, rect_center[1] - 0.5*Ly, 0., rect_center[0] + 0.5*Lx, rect_center[1] + 0.5*Ly + particle_dist + 2*R1, 0.]
inpf.write(" Parameters: " + print_dbl_list(contain_params))

#
# zone info
#
inpf.write("Zone:\n")
inpf.write(" Zones: 2\n")

## zone 1 (bottom particle)
inpf.write(" Zone_1:\n")
inpf.write(" Is_Wall: false\n")

## zone 2 (top particle)
inpf.write(" Zone_2:\n")
inpf.write(" Is_Wall: false\n")

#
# particle info
#
inpf.write("Particle:\n")
inpf.write(" Zone_1:\n")
inpf.write(" Type: rectangle\n")
inpf.write(" Parameters: " + print_dbl_list([Lx, Ly, rect_center[0], rect_center[1], rect_center[2]]))
inpf.write(" Zone_2:\n")
inpf.write(" Type: circle\n")
p2_geom = [R1, cir_center[0], cir_center[1], cir_center[2]]
inpf.write(" Parameters: " + print_dbl_list(p2_geom))

#
# particle generation
#
inpf.write("Particle_Generation:\n")
inpf.write(" From_File: particle_locations_" + str(pp_tag) + ".csv\n")
inpf.write(" File_Data_Type: loc_rad_orient\n")

#
# Mesh info
#
inpf.write("Mesh:\n")

for i in range(len(zones_mesh_fnames)):
inpf.write(" Zone_%d:\n" % (i+1))
inpf.write(" File: %s\n" % (zones_mesh_fnames[i] + "_" + str(pp_tag) + ".msh"))

# Contact info
inpf.write("Contact:\n")

## 11
write_contact_zone_part(inpf, R_contact_factor, damping_active, friction_active, beta_n_eps, friction_coeff, Kn_factor, beta_n_factor, "11", Kn_11)

## copy from 11
copy_contact_zone(inpf, [12, 22], [1, 1])

# Neighbor info
inpf.write("Neighbor:\n")
inpf.write(" Update_Criteria: %s\n" % (neigh_search_criteria))
inpf.write(" Search_Factor: %4.e\n" % (neigh_search_factor))
inpf.write(" Search_Interval: %d\n" % (neigh_search_interval))

# Material info
inpf.write("Material:\n")

## zone 1
write_material_zone_part(inpf, "1", horizon, rho1, K1, G1, Gc1)

## zone 2
inpf.write(" Zone_2:\n")
inpf.write(" Copy_Material_Data: 1\n")

#
# Force
#
if gravity_active == True:
inpf.write("Force_BC:\n")
inpf.write(" Gravity: " + print_dbl_list(gravity))

#
# IC
#
inpf.write("IC:\n")
inpf.write(" Constant_Velocity:\n")
inpf.write(" Velocity_Vector: " + print_dbl_list(free_fall_vel))
inpf.write(" Particle_List: [1]\n")

#
# Displacement
#
inpf.write("Displacement_BC:\n")
inpf.write(" Sets: 1\n")

inpf.write(" Set_1:\n")
inpf.write(" Particle_List: [0]\n")
inpf.write(" Direction: [1,2]\n")
inpf.write(" Time_Function:\n")
inpf.write(" Type: constant\n")
inpf.write(" Parameters:\n")
inpf.write(" - 0.0\n")
inpf.write(" Spatial_Function:\n")
inpf.write(" Type: constant\n")
inpf.write(" Zero_Displacement: true\n")

#
# Output info
#
inpf.write("Output:\n")
inpf.write(" Path: ../out/\n")
inpf.write(" Tags:\n")
inpf.write(" - Displacement\n")
inpf.write(" - Velocity\n")
inpf.write(" - Force\n")
inpf.write(" - Force_Density\n")
inpf.write(" - Damage_Z\n")
inpf.write(" - Damage\n")
inpf.write(" - Nodal_Volume\n")
inpf.write(" - Zone_ID\n")
inpf.write(" - Particle_ID\n")
inpf.write(" - Fixity\n")
inpf.write(" - Force_Fixity\n")
inpf.write(" - Contact_Nodes\n")
inpf.write(" - No_Fail_Node\n")
inpf.write(" - Boundary_Node_Flag\n")
inpf.write(" - Theta\n")
inpf.write(" Output_Interval: %d\n" % (dt_out_n))
inpf.write(" Compress_Type: zlib\n")
inpf.write(" Perform_FE_Out: false\n")
if perform_out:
inpf.write(" Perform_Out: true\n")
else:
inpf.write(" Perform_Out: false\n")
inpf.write(" Test_Output_Interval: %d\n" % (dt_out_n))

inpf.write(" Debug: 3\n")
inpf.write(" Tag_PP: %d\n" %(int(pp_tag)))

# close file
inpf.close()

##-------------------------------------------------------##
##-------------------------------------------------------##
if __name__ == "__main__":
inp_dir = './'
pp_tag = 0
if len(sys.argv) > 1:
pp_tag = int(sys.argv[1])

create_input_file(inp_dir, pp_tag)
Loading