Skip to content

Commit

Permalink
Merge pull request #377 from DongWuTUM/shell_hourglass
Browse files Browse the repository at this point in the history
Modify shell hourglass algorithm and add python script for mutiple runs of test_3d_thin_plate.
  • Loading branch information
Xiangyu-Hu authored Jul 24, 2023
2 parents 5bcdedd + 21785dc commit a4d868f
Show file tree
Hide file tree
Showing 7 changed files with 473 additions and 30 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ ShellStressRelaxationFirstHalf::
gaussian_weight_ = three_gaussian_weights_;
}
/** Define the factor of hourglass control algorithm. */
hourglass_control_factor_ = 0.005;
hourglass_control_factor_ = 0.002;
}
//=================================================================================================//
void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -202,22 +202,24 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
{
Vecd e_ij = inner_neighborhood.e_ij_[n];
Real r_ij = inner_neighborhood.r_ij_[n];
Real dim_inv_r_ij = Dimensions / r_ij;
Real weight = inner_neighborhood.W_ij_[n] * inv_W0_;
Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i],
transformation_matrix_[index_i].transpose() * F_[index_i] * transformation_matrix_[index_i],
pos_[index_j],
transformation_matrix_[index_i].transpose() * F_[index_j] * transformation_matrix_[index_i]);
acceleration += hourglass_control_factor_ * weight * G0_ * pos_jump * dim_inv_r_ij *
inner_neighborhood.dW_ijV_j_[n] * thickness_[index_i];
Real limiter_pos = SMIN(2.0 * pos_jump.norm() / r_ij, 1.0);
acceleration += hourglass_control_factor_ * weight * G0_ * pos_jump * Dimensions *
inner_neighborhood.dW_ijV_j_[n] * limiter_pos;

Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_[index_i] - n0_[index_i],
Vecd pseudo_n_variation_i = pseudo_n_[index_i] - n0_[index_i];
Vecd pseudo_n_variation_j = pseudo_n_[index_j] - n0_[index_j];
Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_variation_i,
transformation_matrix_[index_i].transpose() * F_bending_[index_i] * transformation_matrix_[index_i],
pseudo_n_[index_j] - n0_[index_j],
pseudo_n_variation_j,
transformation_matrix_[index_j].transpose() * F_bending_[index_j] * transformation_matrix_[index_j]);
Vecd rotation_jump = getRotationJump(pseudo_n_jump, transformation_matrix_[index_i]);
pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * rotation_jump * dim_inv_r_ij *
inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 3);
Real limiter_pseudo_n = SMIN(2.0 * pseudo_n_jump.norm() / ((pseudo_n_variation_i- pseudo_n_variation_j).norm() + Eps), 1.0);
pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * pseudo_n_jump * Dimensions *
inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 2) * limiter_pseudo_n;
}

acceleration += (global_stress_i + global_stress_[index_j]) * inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n];
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -173,23 +173,6 @@ Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_
return getWENOStateWithStencilPoints(v1, v2, v3, v4);
}
//=================================================================================================//
Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix)
{
Vec2d local_rotation_jump = Vec2d::Zero();
Vec2d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump;
local_rotation_jump[0] = local_pseuodo_n_jump[0];
return transformation_matrix.transpose() * local_rotation_jump;
}
//=================================================================================================//
Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix)
{
Vec3d local_rotation_jump = Vec3d::Zero();
Vec3d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump;
local_rotation_jump[0] = local_pseuodo_n_jump[0];
local_rotation_jump[1] = local_pseuodo_n_jump[1];
return transformation_matrix.transpose() * local_rotation_jump;
}
//=================================================================================================//
Mat2d getCorrectedAlmansiStrain(const Mat2d &current_local_almansi_strain, const Real &nu_)
{
Mat2d corrected_almansi_strain = current_local_almansi_strain;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -76,10 +76,6 @@ Vecd getWENOLeftState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i
Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value,
const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value);

/** get the artificial rotation from the pseudo-normal jump. */
Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix);
Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix);

/** get the corrected Eulerian Almansi strain tensor according to plane stress problem. */
Mat2d getCorrectedAlmansiStrain(const Mat2d &current_local_almansi_strain, const Real &nu_);
Mat3d getCorrectedAlmansiStrain(const Mat3d &current_local_almansi_strain, const Real &nu_);
Expand Down
27 changes: 27 additions & 0 deletions tests/user_examples/test_3d_thin_plate_python/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
### Python, for ${Python_EXECUTABLE}
find_package(Python3 COMPONENTS Interpreter Development REQUIRED)
### Pybind11
find_package(pybind11 CONFIG REQUIRED)

set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir

set(CMAKE_VERBOSE_MAKEFILE on)

STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} )
PROJECT("${CURRENT_FOLDER}")

SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/")
SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input")
SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload")
SET(BUILD_BIND_PATH "${EXECUTABLE_OUTPUT_PATH}/bind")


file(MAKE_DIRECTORY ${BUILD_BIND_PATH})
file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/pybind_tool/
DESTINATION ${BUILD_BIND_PATH})

aux_source_directory(. DIR_SRCS)
pybind11_add_module(${PROJECT_NAME} ${DIR_SRCS})
set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}")
target_link_libraries(${PROJECT_NAME} PRIVATE sphinxsys_3d)
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
#!/usr/bin/env python3
import os
import sys
import platform
import argparse

# add dynamic link library or shared object to python env
# attention: match current python version with the version exposing the cpp code
sys_str = platform.system()
path_1 = os.path.abspath(os.path.join(os.getcwd(), '../..'))
if sys_str == 'Windows':
path_2 = 'lib'
elif sys_str == 'Linux':
path_2 = 'lib'
else:
# depend on the system
path_2 = 'lib'
path = os.path.join(path_1, path_2)
sys.path.append(path)
# change import depending on the project name
import test_3d_thin_plate_python as test_3d


def run_case(value):
parser = argparse.ArgumentParser()
# set case parameters
parser.add_argument("--restart_step", default=0, type=int)
parser.add_argument("--end_time", default=20, type=int)
parser.add_argument("--loading_factor", default=100, type=float)
case = parser.parse_args()
# project = test_3d.thin_plate_from_sph_cpp(case.loading_factor)
project = test_3d.thin_plate_from_sph_cpp(value)
if project.CmakeTest() == 1:
project.RunCase()
else:
print("check path: ", path)

def net_displacement(file_path, output_file):
with open(file_path, 'r') as infile:
lines = infile.readlines()

first_line_values = lines[1].strip().split()
first_value = float(first_line_values[3])

last_line_values = lines[-1].strip().split()
last_time = float(last_line_values[0])
last_value = float(last_line_values[3])
displacement_z = last_value - first_value

displacement_x = float(last_line_values[1]) - float(first_line_values[1])

with open(output_file, 'a') as outfile:
outfile.write("run_time = " + str(last_time) + ", displacement in Z direction = " + str(displacement_z) + '\n')

return last_time, displacement_z


def copy_files(output_folder):
source_files = ['output/PlateObserver_Position.dat', 'output/SPHBody_PlateBody_0000000000.vtp', 'output/SPHBody_PlateBody_0000000001.vtp']
destination_files = [output_folder + 'PlateObserver_Position.dat', output_folder + 'SPHBody_PlateBody_0000000000.vtp', output_folder + 'SPHBody_PlateBody_0000000001.vtp']

for index, source_file in enumerate(source_files):
# Read the content of the source file
with open(source_file, 'rb') as source:
content = source.read()
# Write the content to the destination file
with open(destination_files[index], 'wb') as destination:
destination.write(content)

def rename_files(value, output_folder):
old_file_name = [output_folder + 'PlateObserver_Position.dat', output_folder + 'SPHBody_PlateBody_0000000000.vtp', output_folder + 'SPHBody_PlateBody_0000000001.vtp']
new_file_name = [(output_folder + 'PlateObserver_Position_' + value + '.dat'), (output_folder + 'SPHBody_PlateBody_0000000000_' + value + '.vtp'), (output_folder + 'SPHBody_PlateBody_0000000001_' + value + '.vtp')]

for index, file in enumerate(old_file_name):
try:
os.rename(file, new_file_name[index])
print('File renamed successfully.')
except OSError:
print('Error: File rename operation failed.')


if __name__ == "__main__":

file_path = 'output/PlateObserver_Position.dat'
output_folder = 'multiple_runs_output/'
output_file = output_folder + 'Displacements.dat'
#values = [6.5, 12.5, 25, 50, 75, 100, 125, 150, 175, 200]
values = [6.5, 12.5]
displacement = list(range(len(values)))
run_time = list(range(len(values)))

if not os.path.exists(output_folder):
# Create the folder
os.makedirs(output_folder)
else:
print(f"The folder '{output_folder}' already exists. Skipping creation.")

with open(output_file, 'w') as outfile:
outfile.truncate()

for index,value in enumerate(values):
with open(output_file, 'a') as outfile:
outfile.write("loading_factor = " + str(value) + '\n')
run_case(value)
run_time[index], displacement[index] = net_displacement(file_path, output_file)
with open(output_file, 'a') as outfile:
outfile.write('\n')
copy_files(output_folder)
rename_files(str(index), output_folder)

print("------------------------------------------------------")
print(f"All the cases finished! Files are saved as follows: ")
for index,value in enumerate(values):
print(f"For loading_factor = {value} :")
new_file_name = [(output_folder + 'PlateObserver_Position_' + str(index) + '.dat'), (output_folder + 'SPHBody_PlateBody_0000000000_' + str(index) + '.vtp'), (output_folder + 'SPHBody_PlateBody_0000000001_' + str(index) + '.vtp')]
print("\t" + new_file_name[0])
print("\t" + new_file_name[1])
print("\t" + new_file_name[2])

print(f"Summary of displacements saved in : " + output_file)
print("\n")
Loading

0 comments on commit a4d868f

Please sign in to comment.