still some issues with this approach

main.py

@@ -11,7 +11,7 @@
 # import projections
 
 # Setup grid parameters
-grid_size = (0.6, 0.2)  # Adjust as needed
+grid_size = (0.15, 0.2)  # Adjust as needed
 cell_size = 0.01  # Adjust based on your simulation scale
 
 # Initialize the grid
@@ -31,14 +31,14 @@
 
 
 # Define disk parameters
-disk1_radius = 8*cell_size
-disk1_center = (0.2, 0.1)
+disk1_radius = 2*cell_size
+disk1_center = (0.05, 0.1)
 disk1_object_id = 1
 
-disk2_radius = 8*cell_size
+disk2_radius = 0*cell_size
 # Calculate the center of the second disk
 # It should be 2 cells from the edge of the first disk
-disk2_x = disk1_center[0] + disk1_radius + 2*cell_size + disk2_radius
+disk2_x = disk1_center[0] + disk1_radius + cell_size + disk2_radius
 disk2_center = (disk2_x, disk1_center[1])
 disk2_object_id = 2
 
@@ -56,13 +56,13 @@
     1: {  # Object ID 1
         "type": "LinearElastic",
         "density": 1000,  # kg/m^3
-        "youngs_modulus": 1e8,  # Pa
+        "youngs_modulus": 1e7,  # Pa
         "poisson_ratio": 0.3
     },
     2: {  # Object ID 2
         "type": "LinearElastic",
         "density": 1000,  # kg/m^3
-        "youngs_modulus": 1e8,  # Pa
+        "youngs_modulus": 1e7,  # Pa
         "poisson_ratio": 0.3
     }
 }
@@ -71,8 +71,8 @@
 particles = Particles(combined_particles, material_properties)
 
 # Set initial velocities for each object
-particles.set_object_velocity(object_id=1, velocity=[0.50, 0.0])  # Object 1 moving right
-particles.set_object_velocity(object_id=2, velocity=[-0.50, 0.0])  # Object 2 moving left
+particles.set_object_velocity(object_id=1, velocity=[1.0, 0.0])  # Object 1 moving right
+particles.set_object_velocity(object_id=2, velocity=[-1.0, 0.0])  # Object 2 moving left
 
 print(f"Generated {len(particles1)} particles for disk 1")
 print(f"Generated {len(particles2)} particles for disk 2")
@@ -93,7 +93,7 @@
     else:
         print("Invalid input. Please press Enter to continue or '0' to exit.")
 
-dt = 1e-6#particles.compute_dt()
+dt = 1e-9#particles.compute_dt()
 
 print(f"dt: {dt}")
 # Create MPM solver
@@ -102,20 +102,20 @@
 
 
 
-num_steps = 10000  # TODO: Adjust as needed
+num_steps = 100000  # TODO: Adjust as needed
 
 
 
 # Create an output directory if it doesn't exist
-output_dir = "simulation_output"
+output_dir = "simulation_output_low_res"
 os.makedirs(output_dir, exist_ok=True)
 
 # Main simulation loop
 for step in tqdm(range(num_steps), desc="Simulating"):
     solver.step()
 
     # Save outputs every 100 steps
-    if step % 100 == 0:
+    if step % 1000 == 0:
         output_file = os.path.join(output_dir, f"step_{step:05d}.csv")
         save_particles_to_csv(particles, output_file)
         print(f"Saved output for step {step} to {output_file}")

material.py

@@ -5,23 +5,25 @@ def __init__(self, density, youngs_modulus, poisson_ratio):
         self.density = density
         self.youngs_modulus = youngs_modulus
         self.poisson_ratio = poisson_ratio
-
-        # Compute Lame parameters
-        self.lambda_ = (youngs_modulus * poisson_ratio) / ((1 + poisson_ratio) * (1 - 2 * poisson_ratio))
-        self.mu = youngs_modulus / (2 * (1 + poisson_ratio))
 
     def compute_stress_rate(self, strain_rate):
-        # Compute stress rate using Hooke's law for 2D stress state
-        # Expect strain_rate as a vector [exx, eyy, exy]
         assert strain_rate.shape == (2,2), "Strain rate must be a 2x2 matrix"
 
-        exx, eyy, exy, eyx = strain_rate.flatten()
+        # Calculate bulk modulus K and shear modulus G from Young's modulus and Poisson's ratio
+        K = self.youngs_modulus / (3 * (1 - 2 * self.poisson_ratio))
+        G = self.youngs_modulus / (2 * (1 + self.poisson_ratio))
+
+        # Calculate volumetric strain rate
+        volumetric_strain_rate = np.trace(strain_rate)
+
+        # Calculate deviatoric strain rate
+        deviatoric_strain_rate = strain_rate - (volumetric_strain_rate / 3) * np.eye(2)
+
+        # Compute stress rate using volumetric and deviatoric parts
+        volumetric_stress_rate = 3 * K * volumetric_strain_rate * np.eye(2)
+        deviatoric_stress_rate = 2 * G * deviatoric_strain_rate
 
-        # Compute stress components
-        factor = self.youngs_modulus / ((1 + self.poisson_ratio) * (1 - 2 * self.poisson_ratio))
-        sxx = factor * ((1 - self.poisson_ratio) * exx + self.poisson_ratio * (eyy))
-        syy = factor * ((1 - self.poisson_ratio) * eyy + self.poisson_ratio * (exx))
-        sxy = self.mu * 2 * exy
+        stress_rate = volumetric_stress_rate + deviatoric_stress_rate
 
-        return np.array([[sxx, sxy], [sxy, syy]])
+        return stress_rate
 

nodes.py

@@ -52,6 +52,7 @@ def __init__(self, position):
 
     def reset(self):
         self.mass = 0.0
+        self.momentum.fill(0)
         self.velocity.fill(0)
         self.force.fill(0)
         # Reset other properties as needed

particle_shapes.py

@@ -15,9 +15,9 @@ def _generate_grid_particles(self, x_range, y_range):
             for j in range(y_range[0], y_range[1]):
                 for px in range(2):
                     for py in range(2):
-                        # Place particles at 1/4 and 3/4 of the cell in each dimension
-                        x = (i + (px + 1) / 3) * self.cell_size
-                        y = (j + (py + 1) / 3) * self.cell_size
+                        # Place particles at the center of each quadrant (1/4 and 3/4 of the cell)
+                        x = (i + (px + 0.5) / 2) * self.cell_size
+                        y = (j + (py + 0.5) / 2) * self.cell_size
 
                         if self._is_inside(x, y):
                             particle = {
@@ -70,8 +70,9 @@ def generate_particles(self):
             for j in range(int(y_min / self.cell_size), int(y_max / self.cell_size) + 1):
                 for px in range(2):
                     for py in range(2):
-                        x = (i + (px + 1) / 3) * self.cell_size
-                        y = (j + (py + 1) / 3) * self.cell_size
+                        # Place particles at the center of each quadrant (1/4 and 3/4 of the cell)
+                        x = (i + (px + 0.5) / 2) * self.cell_size
+                        y = (j + (py + 0.5) / 2) * self.cell_size
                         if self._is_inside(x, y):
                             particle = {
                                 'position': np.array([x, y]),

particles.py

@@ -27,11 +27,13 @@ def create_particles(self, particle_data):
                 particle = {
                     'position': position,
                     'velocity': velocity,
+                    'Gvelocity': np.zeros(2),
                     'acceleration': np.zeros(2),
                     'strain_rate': np.zeros((2, 2)),
                     'density': material_props['density'],
                     'mass': data['mass'],
-                    'volume': data['volume'],
+                    'volume': data['mass']/material_props['density'],
+                    'density_rate': np.zeros(2),
                     'stress': np.zeros((2, 2)),
                     'strain': np.zeros((2, 2)),
                     'id': data.get('id', len(self.particles)),

results_vis.py

@@ -4,12 +4,15 @@
 import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.colors import Normalize
+import argparse
 
-# Directory containing the CSV files
-data_dir = "simulation_output"
+# Add command-line argument parsing
+parser = argparse.ArgumentParser(description="Visualize simulation results.")
+parser.add_argument("data_dir", help="Directory containing the CSV files")
+args = parser.parse_args()
 
-# Get all CSV files in the directory
-csv_files = sorted([f for f in os.listdir(data_dir) if f.endswith('.csv')])
+# Use the provided data directory
+data_dir = args.data_dir
 
 # Variables to store grid boundaries
 grid_min_x = float('inf')
@@ -42,6 +45,7 @@ def read_csv(filename):
     return positions, velocities, stress, volumes, masses
 
 # Read the first file to get the number of particles and initialize the plot
+csv_files = sorted([f for f in os.listdir(data_dir) if f.endswith('.csv')])
 initial_positions, initial_velocities, initial_stress, initial_volumes, initial_masses = read_csv(csv_files[0])
 num_particles = len(initial_positions)
 
@@ -52,8 +56,8 @@ def read_csv(filename):
 ax2 = fig.add_subplot(gs[0, 1])
 ax3 = fig.add_subplot(gs[1, :])
 
-scatter1 = ax1.scatter([], [], s=10, c=[], cmap='viridis')
-scatter2 = ax2.scatter([], [], s=10, c=[], cmap='plasma')
+scatter1 = ax1.scatter([], [], s=50, c=[], cmap='viridis', marker='s')
+scatter2 = ax2.scatter([], [], s=50, c=[], cmap='plasma', marker='s')
 line1, = ax3.plot([], [], label='Kinetic Energy')
 line2, = ax3.plot([], [], label='Elastic Energy')
 line3, = ax3.plot([], [], label='Total Energy')
@@ -129,13 +133,13 @@ def update(frame):
     stress_max = max(stress_max, von_mises.max())
 
     # Update color normalizations
-    # norm1.autoscale(vel_mag)
+    norm1.autoscale(vel_mag)
     # norm1.autoscale(velocities[:, 0])
-    norm1.vmin = -2.0#vel_mag.min()
-    norm1.vmax = 2.0
-    # norm2.autoscale(von_mises)
-    norm2.vmin = 0.0#stress_min
-    norm2.vmax = 1.0e5#stress_max
+    # norm1.vmin = -2.0#vel_mag.min()
+    # norm1.vmax = 2.0
+    norm2.autoscale(von_mises)
+    # norm2.vmin = 0.0#stress_min
+    # norm2.vmax = 1.0e5#stress_max
 
     # scatter1.set_array(vel_mag)
     scatter1.set_array(velocities[:, 0])

shape_function.py

@@ -1,4 +1,5 @@
 import numpy as np
+import matplotlib.pyplot as plt
 
 def shape_function(particle_pos, node_pos, cell_size):
     """
@@ -73,4 +74,34 @@ def piecewise_gradient(xp, xv):
     dsy_dy = piecewise_gradient(particle_pos[1], node_pos[1])
 
     return dsx_dx, dsy_dy
+# wrap in main
+if __name__ == "__main__":  
+    # Parameters
+    cell_size = 1.0
+    node_pos = np.array([0.0, 0.0])
+    particle_positions = np.linspace(-1.5 * cell_size, 1.5 * cell_size, 500)
+
+    # Calculate shape function and gradient values
+    shape_values = [shape_function(np.array([x, 0.0]), node_pos, cell_size)[0] for x in particle_positions]
+    gradient_values = [gradient_shape_function(np.array([x, 0.0]), node_pos, cell_size)[0] for x in particle_positions]
+
+    # Plotting
+    plt.figure(figsize=(12, 6))
+
+    plt.subplot(1, 2, 1)
+    plt.plot(particle_positions, shape_values, label='Shape Function')
+    plt.title('Shape Function')
+    plt.xlabel('Particle Position')
+    plt.ylabel('Shape Function Value')
+    plt.legend()
+
+    plt.subplot(1, 2, 2)
+    plt.plot(particle_positions, gradient_values, label='Gradient of Shape Function', color='r')
+    plt.title('Gradient of Shape Function')
+    plt.xlabel('Particle Position')
+    plt.ylabel('Gradient Value')
+    plt.legend()
+
+    plt.tight_layout()
+    plt.show()