Modify shell hourglass algorithm and add python script for mutiple runs of test_3d_thin_plate.

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp

@@ -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)

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h

@@ -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];

src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp

@@ -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;

src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h

@@ -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_);

tests/user_examples/test_3d_thin_plate_python/CMakeLists.txt

@@ -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)

tests/user_examples/test_3d_thin_plate_python/pybind_tool/pybind_test.py

@@ -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")