Merge pull request #1 from SamRaymond/c-like

C like

README.md

@@ -41,6 +41,10 @@ The Material Point Method is a numerical technique used to simulate the behavior
 8. **Results Processing (`results.py`)**: 
    - Handles saving simulation data to CSV files
 
+9. **Particle Shapes (`particle_shapes.py`)**: 
+   - Defines different particle shapes (e.g., Disk, Block)
+   - Generates particles for each shape
+
 ## Key Features
 
 - Simulation of two colliding elastic disks
@@ -49,6 +53,7 @@ The Material Point Method is a numerical technique used to simulate the behavior
 - Background grid for computation
 - Visualization of particle positions, velocities, and stresses
 - Energy tracking (kinetic, elastic, and total)
+- Flexible particle shape generation (currently supports Disk and Block)
 
 ## Usage
 
@@ -62,6 +67,7 @@ The Material Point Method is a numerical technique used to simulate the behavior
 - Modify grid size and resolution in `main.py`
 - Implement new material models in `material.py`
 - Add boundary conditions in `solver.py`
+- Create new particle shapes in `particle_shapes.py`
 
 ## Output
 
@@ -88,6 +94,7 @@ The `results_vis.py` script provides an animated visualization of:
 - The current implementation focuses on 2D simulations
 - Only linear elastic materials are implemented, but the structure allows for easy addition of other material models
 - The simulation uses an explicit time integration scheme
+- Particle shapes can be easily extended by subclassing the `ParticleShape` class
 
 ## Future Improvements
 
@@ -96,6 +103,7 @@ The `results_vis.py` script provides an animated visualization of:
 - Implement adaptive time-stepping
 - Optimize performance for larger simulations
 - Add more boundary condition options
+- Implement additional particle shapes
 
 ## Dependencies
 
@@ -107,4 +115,8 @@ The `results_vis.py` script provides an animated visualization of:
 
 1. Ensure all dependencies are installed
 2. Run `python main.py` to start the simulation
-3. After the simulation completes, run `python results_vis.py` to visualize the results
+3. After the simulation completes, run `python results_vis.py <output_directory>` to visualize the results
+
+## Visualization Options
+
+The `results_vis.py` script now accepts command-line arguments for specifying the data directory. Use it as follows:

main.py

@@ -11,7 +11,7 @@
 # import projections
 
 # Setup grid parameters
-grid_size = (0.15, 0.2)  # Adjust as needed
+grid_size = (0.5, 0.5)  # Adjust as needed
 cell_size = 0.01  # Adjust based on your simulation scale
 
 # Initialize the grid
@@ -27,97 +27,96 @@
 # For example, to access the node at position (i, j):
 # node = grid.get_node(i, j)
 
-# Initialize Material Points
 
+object_ids = [1,2]
+
+# Initialize Material Points
+# Define material properties
+material_properties = {
+    object_ids[0]: {  # Object ID 1
+        "type": "LinearElastic",
+        "density": 1000,  # kg/m^3
+        "youngs_modulus": 1e8,  # Pa
+        "poisson_ratio": 0.3
+    },
+    object_ids[1]: {  # Object ID 2
+        "type": "LinearElastic",
+        "density": 1000,  # kg/m^3
+        "youngs_modulus": 1e8,  # Pa
+        "poisson_ratio": 0.3
+    }
+}
 
 # Define disk parameters
-disk1_radius = 2*cell_size
-disk1_center = (0.05, 0.1)
+disk1_radius = 6*cell_size
+disk1_center = (0.25, 0.25)
 disk1_object_id = 1
 
-disk2_radius = 0*cell_size
+disk2_radius = 6*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 + cell_size + disk2_radius
+disk2_x = disk1_center[0] + disk1_radius + 3*cell_size + disk2_radius
 disk2_center = (disk2_x, disk1_center[1])
 disk2_object_id = 2
 
 # Create Disk objects
-disk1 = Disk(cell_size, disk1_radius, disk1_center, disk1_object_id)
-disk2 = Disk(cell_size, disk2_radius, disk2_center, disk2_object_id)
+disk1 = Disk(cell_size, disk1_radius, disk1_center, object_ids[0], material_properties[object_ids[0]]["density"])
+disk2 = Disk(cell_size, disk2_radius, disk2_center, object_ids[1], material_properties[object_ids[1]]["density"])
 
 # Generate particles for both disks
 particles1 = disk1.generate_particles()
 particles2 = disk2.generate_particles()
 combined_particles = particles1 + particles2
 
-# Define material properties
-material_properties = {
-    1: {  # Object ID 1
-        "type": "LinearElastic",
-        "density": 1000,  # kg/m^3
-        "youngs_modulus": 1e7,  # Pa
-        "poisson_ratio": 0.3
-    },
-    2: {  # Object ID 2
-        "type": "LinearElastic",
-        "density": 1000,  # kg/m^3
-        "youngs_modulus": 1e7,  # Pa
-        "poisson_ratio": 0.3
-    }
-}
+
 
 # Create Particles object
-particles = Particles(combined_particles, material_properties)
+particles = Particles(combined_particles, material_properties,cell_size)
 
 # Set initial velocities for each object
-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
+particles.set_object_velocity(object_id=object_ids[0], velocity=[1.0, 0.0])  # Object 1 moving right
+particles.set_object_velocity(object_id=object_ids[1], 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")
 print(f"Total particles: {particles.get_particle_count()}")
 
-# # Call the visualization function
-visualize_particles_and_grid(particles, grid, 0)
+# # # Call the visualization function
+# visualize_particles_and_grid(particles, grid, 0)
 
-# Wait for keyboard input to continue or exit
-while True:
-    user_input = input("Press Enter to continue or '0' to exit: ")
-    if user_input == '0':
-        print("Exiting simulation...")
-        exit()
-    elif user_input == '':
-        print("Continuing simulation...")
-        break
-    else:
-        print("Invalid input. Please press Enter to continue or '0' to exit.")
+# # Wait for keyboard input to continue or exit
+# while True:
+#     user_input = input("Press Enter to continue or '0' to exit: ")
+#     if user_input == '0':
+#         print("Exiting simulation...")
+#         exit()
+#     elif user_input == '':
+#         print("Continuing simulation...")
+#         break
+#     else:
+#         print("Invalid input. Please press Enter to continue or '0' to exit.")
 
-dt = 1e-9#particles.compute_dt()
+dt = particles.compute_dt(cfl_factor=0.2)
 
 print(f"dt: {dt}")
 # Create MPM solver
 
 solver = MPMSolver(particles, grid, dt)
 
-
-
-num_steps = 100000  # TODO: Adjust as needed
-
-
+num_steps = 10000
 
 # Create an output directory if it doesn't exist
-output_dir = "simulation_output_low_res"
+output_dir = "simulation_output"
 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 % 1000 == 0:
+    if step % 100 == 0:
         output_file = os.path.join(output_dir, f"step_{step:05d}.csv")
-        save_particles_to_csv(particles, output_file)
+        save_particles_to_csv(particles,output_file)
         print(f"Saved output for step {step} to {output_file}")
 
 print("Simulation completed successfully!")