iniital upload - working collision of two disks

README.md

@@ -1 +1,110 @@
-# py_mpm
+# MPM Simulation
+
+This project implements a 2D Material Point Method (MPM) simulation for elastic materials. It simulates the interaction between two circular disks using MPM techniques.
+
+## Overview
+
+The Material Point Method is a numerical technique used to simulate the behavior of solids, fluids, and other materials. This implementation focuses on simulating elastic materials in a 2D space.
+
+## Key Components
+
+1. **Main Simulation (`main.py`)**: 
+   - Sets up the simulation parameters
+   - Initializes the grid and particles
+   - Runs the main simulation loop
+
+2. **MPM Solver (`solver.py`)**: 
+   - Implements the core MPM algorithm
+   - Handles particle-to-grid and grid-to-particle transfers
+   - Updates particle and grid properties each time step
+
+3. **Particles (`particles.py`)**: 
+   - Manages particle data and properties
+   - Handles material assignment to particles
+
+4. **Grid (`nodes.py`)**: 
+   - Defines the background grid structure
+   - Manages grid nodes and their properties
+
+5. **Material Models (`material.py`)**: 
+   - Implements material behavior (currently Linear Elastic)
+   - Computes stress based on strain
+
+6. **Shape Functions (`shape_function.py`)**: 
+   - Defines shape functions for MPM interpolation
+   - Includes functions for shape function gradients
+
+7. **Visualization (`plotting.py`, `results_vis.py`)**: 
+   - Provides functions to visualize the simulation results
+   - Includes both static plotting and animation capabilities
+
+8. **Results Processing (`results.py`)**: 
+   - Handles saving simulation data to CSV files
+
+## Key Features
+
+- Simulation of two colliding elastic disks
+- Linear elastic material model
+- Particle-based representation of materials
+- Background grid for computation
+- Visualization of particle positions, velocities, and stresses
+- Energy tracking (kinetic, elastic, and total)
+
+## Usage
+
+1. Set up the simulation parameters in `main.py`
+2. Run `main.py` to start the simulation
+3. Use `results_vis.py` to visualize the simulation results
+
+## Customization
+
+- Adjust material properties in `main.py`
+- Modify grid size and resolution in `main.py`
+- Implement new material models in `material.py`
+- Add boundary conditions in `solver.py`
+
+## Output
+
+The simulation generates CSV files for each time step, storing particle data including:
+- Position
+- Velocity
+- Stress
+- Strain
+- Volume
+- Mass
+- Density
+- Material properties
+
+## Visualization
+
+The `results_vis.py` script provides an animated visualization of:
+- Particle positions
+- Velocity magnitudes
+- Von Mises stress
+- Energy evolution over time
+
+## Notes
+
+- 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
+
+## Future Improvements
+
+- Implement more advanced material models (e.g., plasticity, damage)
+- Add support for 3D simulations
+- Implement adaptive time-stepping
+- Optimize performance for larger simulations
+- Add more boundary condition options
+
+## Dependencies
+
+- NumPy
+- Matplotlib
+- tqdm (for progress bars)
+
+## Running the Simulation
+
+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

main.py

@@ -0,0 +1,123 @@
+import os
+import csv
+import numpy as np
+from tqdm import tqdm
+from particle_shapes import Disk
+from nodes import Grid
+from particles import Particles
+from solver import MPMSolver
+from plotting import visualize_particles_and_grid
+from results import save_particles_to_csv
+# import projections
+
+# Setup grid parameters
+grid_size = (0.6, 1)  # Adjust as needed
+cell_size = 0.01  # Adjust based on your simulation scale
+
+# Initialize the grid
+grid = Grid(grid_size, cell_size)
+
+# Calculate the physical size of the grid
+physical_size = (grid_size[0] * cell_size, grid_size[1] * cell_size)
+
+print(f"Grid initialized with size {grid_size} and physical dimensions {physical_size}")
+print(f"Total number of nodes: {grid.total_nodes}")  # We now have 2500 nodes (50 * 50)
+
+# Nodes can be accessed using grid.nodes
+# For example, to access the node at position (i, j):
+# node = grid.get_node(i, j)
+
+# Initialize Material Points
+
+
+# Define disk parameters
+disk1_radius = 8*cell_size
+disk1_center = (0.2, 0.5)
+disk1_object_id = 1
+
+disk2_radius = 8*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 + 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)
+
+# 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": 1e6,  # Pa
+        "poisson_ratio": 0.3
+    },
+    2: {  # Object ID 2
+        "type": "LinearElastic",
+        "density": 1000,  # kg/m^3
+        "youngs_modulus": 1e6,  # Pa
+        "poisson_ratio": 0.3
+    }
+}
+
+# Create Particles object
+particles = Particles(combined_particles, material_properties)
+
+# Set initial velocities for each object
+particles.set_object_velocity(object_id=1, velocity=[2.0, 0.0])  # Object 1 moving right
+particles.set_object_velocity(object_id=2, velocity=[-2.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)
+
+# 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-6#particles.compute_dt()
+
+print(f"dt: {dt}")
+# Create MPM solver
+
+solver = MPMSolver(particles, grid, dt)
+
+
+
+num_steps = 10000  # TODO: Adjust as needed
+
+
+
+# Create an output directory if it doesn't exist
+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:
+        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}")
+
+print("Simulation completed successfully!")

material.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+class LinearElastic:
+    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()
+
+        # 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
+
+        return np.array([[sxx, sxy], [sxy, syy]])
+

nodes.py

@@ -0,0 +1,57 @@
+import numpy as np
+
+class Grid:
+    def __init__(self, physical_size, cell_size):
+        self.physical_size = np.array(physical_size)
+        self.cell_size = cell_size
+        self.node_count = np.ceil(self.physical_size / self.cell_size).astype(int) + 1
+        self.nodes = self.initialize_grid()
+        self.total_nodes = np.prod(self.node_count)
+
+    def initialize_grid(self):
+        grid = np.empty(self.node_count, dtype=object)
+        for i in range(self.node_count[0]):
+            for j in range(self.node_count[1]):
+                x = i * self.cell_size
+                y = j * self.cell_size
+                grid[i, j] = Node(position=np.array([x, y]))
+        return grid
+
+    def reset_nodes(self):
+        for i in range(self.node_count[0]):
+            for j in range(self.node_count[1]):
+                self.nodes[i, j].reset()
+
+    def get_node(self, i, j):
+        return self.nodes[i, j]
+
+    def set_node(self, i, j, node):
+        self.nodes[i, j] = node
+
+    def get_nearby_nodes(self, particle_position):
+        i = int(particle_position[0] / self.cell_size)
+        j = int(particle_position[1] / self.cell_size)
+        nearby_nodes = []
+        for di in range(-1, 2):
+            for dj in range(-1, 2):
+                if 0 <= i + di < self.node_count[0] and 0 <= j + dj < self.node_count[1]:
+                    nearby_nodes.append(self.nodes[i + di, j + dj])
+        return nearby_nodes
+
+    @property
+    def node_count_total(self):
+        return self.total_nodes
+
+class Node:
+    def __init__(self, position):
+        self.position = position
+        self.mass = 0.0
+        self.velocity = np.zeros(2)
+        self.momentum = np.zeros(2)
+        self.force = np.zeros(2)
+
+    def reset(self):
+        self.mass = 0.0
+        self.velocity.fill(0)
+        self.force.fill(0)
+        # Reset other properties as needed