I will implement a Laplace solver for scattering off two touching cir…

geometry.py

@@ -0,0 +1,73 @@
+import numpy as np
+
+# Define radii
+r1 = 1.0
+r2 = 1.0
+
+# Define centers
+center1 = (-r1, 0)
+center2 = (r2, 0)
+
+def parameterize_circle(center, radius, num_points):
+    """
+    Generates points on the circumference of a circle.
+
+    Args:
+        center (tuple): (x, y) coordinates of the circle's center.
+        radius (float): The radius of the circle.
+        num_points (int): The number of points to generate.
+
+    Returns:
+        tuple: Two NumPy arrays (x_coords, y_coords) representing the points.
+    """
+    angles = np.linspace(0, 2 * np.pi, num_points, endpoint=False)
+    x_coords = center[0] + radius * np.cos(angles)
+    y_coords = center[1] + radius * np.sin(angles)
+    return x_coords, y_coords
+
+def incident_field(x, y):
+    """
+    Defines the incident field.
+    For now, phi_incident(x, y) = x.
+    """
+    return x
+
+def boundary_condition_on_surface(x, y, incident_field_func):
+    """
+    Calculates the desired value for the scattered field on the surface.
+    This is -incident_field_func(x, y).
+    """
+    return -incident_field_func(x, y)
+
+if __name__ == "__main__":
+    num_points_c1 = 100
+    num_points_c2 = 100
+
+    x_c1, y_c1 = parameterize_circle(center1, r1, num_points_c1)
+    x_c2, y_c2 = parameterize_circle(center2, r2, num_points_c2)
+
+    print("First 5 points of C1:")
+    for i in range(5):
+        print(f"({x_c1[i]:.2f}, {y_c1[i]:.2f})")
+
+    print("\nFirst 5 points of C2:")
+    for i in range(5):
+        print(f"({x_c2[i]:.2f}, {y_c2[i]:.2f})")
+
+    # Calculate boundary conditions for C1
+    bc_c1 = np.zeros(num_points_c1)
+    for i in range(num_points_c1):
+        bc_c1[i] = boundary_condition_on_surface(x_c1[i], y_c1[i], incident_field)
+
+    print("\nFirst 5 boundary condition values for C1:")
+    for i in range(5):
+        print(f"{bc_c1[i]:.2f}")
+
+    # Calculate boundary conditions for C2
+    bc_c2 = np.zeros(num_points_c2)
+    for i in range(num_points_c2):
+        bc_c2[i] = boundary_condition_on_surface(x_c2[i], y_c2[i], incident_field)
+
+    print("\nFirst 5 boundary condition values for C2:")
+    for i in range(5):
+        print(f"{bc_c2[i]:.2f}")

solver.py

@@ -0,0 +1,118 @@
+import numpy as np
+from geometry import parameterize_circle, incident_field, boundary_condition_on_surface
+from solver_setup import place_poles, construct_system_matrix
+
+# Define problem parameters (can be adjusted)
+r1_val = 1.0
+r2_val = 1.0
+center1_val = (-r1_val, 0)
+center2_val = (r2_val, 0)
+
+num_poles_c = 50  # Number of poles per circle
+pole_sf = 0.95    # Pole scale factor
+num_coll_pts_c = 50 # Number of collocation points per circle
+
+def setup_problem(r1, r2, center1, center2, num_poles_circle, pole_scale_factor, num_collocation_points_circle):
+    """
+    Generates all necessary components for the scattering problem.
+    """
+    # Generate Collocation Points for C1 and C2
+    coll_x1, coll_y1 = parameterize_circle(center1, r1, num_collocation_points_circle)
+    coll_x2, coll_y2 = parameterize_circle(center2, r2, num_collocation_points_circle)
+    all_coll_x = np.concatenate((coll_x1, coll_x2))
+    all_coll_y = np.concatenate((coll_y1, coll_y2))
+    z_collocation_pts = all_coll_x + 1j * all_coll_y
+
+    # Generate Pole Locations for C1 and C2
+    poles_x1, poles_y1 = place_poles(center1, r1, num_poles_circle, pole_scale_factor)
+    poles_x2, poles_y2 = place_poles(center2, r2, num_poles_circle, pole_scale_factor)
+    all_poles_x = np.concatenate((poles_x1, poles_x2))
+    all_poles_y = np.concatenate((poles_y1, poles_y2))
+    z_poles_pts = all_poles_x + 1j * all_poles_y
+
+    # Calculate Boundary Values (b_vector) at Collocation Points
+    b_vec = np.zeros(len(all_coll_x))
+    for i in range(len(all_coll_x)):
+        # Using boundary_condition_on_surface which is -incident_field
+        b_vec[i] = boundary_condition_on_surface(all_coll_x[i], all_coll_y[i], incident_field)
+
+    # Construct System Matrix A
+    A_mat = construct_system_matrix(z_collocation_pts, z_poles_pts)
+
+    return z_collocation_pts, z_poles_pts, A_mat, b_vec
+
+# Solve for Pole Coefficients
+# This will be done in the main block after calling setup_problem
+
+def calculate_scattered_field(x, y, z_p, coeffs):
+    """
+    Calculates the scattered field at a point (x,y) due to poles.
+    """
+    z_eval = x + 1j * y
+    scattered_field_value = 0.0
+    for k in range(len(z_p)):
+        if not np.isclose(z_eval, z_p[k]): # Avoid division by zero if poles are at eval point
+            scattered_field_value += np.real(coeffs[k] / (z_eval - z_p[k]))
+    return scattered_field_value
+
+def calculate_total_field(x, y, z_p, coeffs, inc_field_func):
+    """
+    Calculates the total field (incident + scattered) at a point (x,y).
+    """
+    scattered_val = calculate_scattered_field(x, y, z_p, coeffs)
+    incident_val = inc_field_func(x, y)
+    return scattered_val + incident_val
+
+if __name__ == "__main__":
+    print("Setting up the scattering problem...")
+    z_collocation, z_poles, A_matrix, b_vector = setup_problem(
+        r1_val, r2_val, center1_val, center2_val,
+        num_poles_c, pole_sf, num_coll_pts_c
+    )
+
+    print(f"A_matrix shape: {A_matrix.shape}")
+    print(f"b_vector shape: {b_vector.shape}")
+    print(f"z_poles shape: {z_poles.shape}")
+    print(f"z_collocation shape: {z_collocation.shape}")
+
+    print("\nSolving for pole coefficients...")
+    coefficients = np.linalg.solve(A_matrix, b_vector)
+
+    print("\nFirst 10 coefficients:")
+    for i in range(min(10, len(coefficients))):
+        print(f"Coeff {i}: {coefficients[i]:.4f}")
+
+    test_points = [
+        (2.0, 0.0),
+        (0.0, 2.0),
+        (-2.0, 0.5), # Test point near C1 but outside
+        (2.0, -0.5)  # Test point near C2 but outside
+    ]
+
+    print("\nField calculations at test points (outside circles):")
+    for pt_x, pt_y in test_points:
+        sc_field = calculate_scattered_field(pt_x, pt_y, z_poles, coefficients)
+        tot_field = calculate_total_field(pt_x, pt_y, z_poles, coefficients, incident_field)
+        print(f"Point ({pt_x:.1f}, {pt_y:.1f}): Scattered = {sc_field:.4f}, Total = {tot_field:.4f}, Incident = {incident_field(pt_x, pt_y):.4f}")
+
+    print("\nVerification: Total field at first 3 collocation points (should be ~0):")
+    for i in range(min(3, len(z_collocation))):
+        coll_x = np.real(z_collocation[i])
+        coll_y = np.imag(z_collocation[i])
+        tot_field_coll = calculate_total_field(coll_x, coll_y, z_poles, coefficients, incident_field)
+        # Expected b_vector value at this point (which is -incident_field)
+        expected_boundary_val = b_vector[i]
+        # Total field = scattered + incident. Boundary condition was scattered = -incident. So total = 0.
+        print(f"Collocation Pt ({coll_x:.2f}, {coll_y:.2f}): Total Field = {tot_field_coll:.4e} (Expected on boundary: ~0, b_vector value: {expected_boundary_val:.2f})")
+
+    print("\nVerification: Total field at a point on C1 (idx N_coll_c // 2) and C2 (idx N_coll_c + N_coll_c // 2) (should be ~0):")
+    if len(z_collocation) > num_coll_pts_c + num_coll_pts_c // 2:
+        indices_to_check = [num_coll_pts_c // 2, num_coll_pts_c + num_coll_pts_c // 2 -1] # one from each circle
+        for idx in indices_to_check:
+            coll_x = np.real(z_collocation[idx])
+            coll_y = np.imag(z_collocation[idx])
+            tot_field_coll = calculate_total_field(coll_x, coll_y, z_poles, coefficients, incident_field)
+            expected_boundary_val = b_vector[idx]
+            print(f"Collocation Pt ({coll_x:.2f}, {coll_y:.2f}): Total Field = {tot_field_coll:.4e} (Expected on boundary: ~0, b_vector value: {expected_boundary_val:.2f})")
+    else:
+        print("Not enough collocation points to run the second verification set.")

solver_setup.py

@@ -0,0 +1,113 @@
+import numpy as np
+from geometry import r1, r2, center1, center2, parameterize_circle, incident_field, boundary_condition_on_surface
+
+# Define Pole Placement Parameters
+num_poles_circle = 50  # number of poles per circle
+pole_scale_factor = 0.95  # scales poles inward from the boundary
+
+def place_poles(center, radius, num_poles, scale_factor):
+    """
+    Places poles inside a circle, scaled inward from the boundary.
+
+    Args:
+        center (tuple): (x, y) coordinates of the circle's center.
+        radius (float): The radius of the circle.
+        num_poles (int): The number of poles to generate.
+        scale_factor (float): Factor to scale poles inward from the boundary.
+
+    Returns:
+        tuple: Two NumPy arrays (pole_x_coords, pole_y_coords).
+    """
+    # Generate points on the boundary
+    boundary_x, boundary_y = parameterize_circle(center, radius, num_poles)
+
+    # Scale points inward to get pole locations
+    pole_x_coords = center[0] + (boundary_x - center[0]) * scale_factor
+    pole_y_coords = center[1] + (boundary_y - center[1]) * scale_factor
+
+    return pole_x_coords, pole_y_coords
+
+# Generate Collocation Points
+num_collocation_points_circle = 50  # can be same as num_poles for now
+coll_x1, coll_y1 = parameterize_circle(center1, r1, num_collocation_points_circle)
+coll_x2, coll_y2 = parameterize_circle(center2, r2, num_collocation_points_circle)
+
+all_coll_x = np.concatenate((coll_x1, coll_x2))
+all_coll_y = np.concatenate((coll_y1, coll_y2))
+z_collocation = all_coll_x + 1j * all_coll_y
+
+# Generate Pole Locations
+poles_x1, poles_y1 = place_poles(center1, r1, num_poles_circle, pole_scale_factor)
+poles_x2, poles_y2 = place_poles(center2, r2, num_poles_circle, pole_scale_factor)
+
+all_poles_x = np.concatenate((poles_x1, poles_x2))
+all_poles_y = np.concatenate((poles_y1, poles_y2)) # Corrected concatenation
+z_poles = all_poles_x + 1j * all_poles_y
+
+# Define System Matrix Construction Function
+def construct_system_matrix(z_coll, z_p):
+    """
+    Constructs the system matrix A for the method of fundamental solutions.
+    A[i, j] = real(1 / (z_coll[i] - z_p[j]))
+
+    Args:
+        z_coll (np.array): Complex array of collocation points.
+        z_p (np.array): Complex array of pole locations.
+
+    Returns:
+        np.array: The system matrix A.
+    """
+    num_coll_points = len(z_coll)
+    num_poles = len(z_p)
+    A_matrix = np.zeros((num_coll_points, num_poles))
+
+    for i in range(num_coll_points):
+        for j in range(num_poles):
+            A_matrix[i, j] = np.real(1 / (z_coll[i] - z_p[j]))
+    return A_matrix
+
+# Calculate Boundary Values at Collocation Points
+b_vector = np.zeros(len(all_coll_x))
+for i in range(len(all_coll_x)):
+    incident_val = incident_field(all_coll_x[i], all_coll_y[i])
+    b_vector[i] = -incident_val # Using the simplified direct relation for now
+
+if __name__ == "__main__":
+    print("Shapes of generated arrays:")
+    print(f"all_coll_x: {all_coll_x.shape}")
+    print(f"z_collocation: {z_collocation.shape}")
+    print(f"all_poles_x: {all_poles_x.shape}")
+    print(f"z_poles: {z_poles.shape}")
+    print(f"b_vector: {b_vector.shape}")
+
+    # Construct System Matrix
+    A_matrix = construct_system_matrix(z_collocation, z_poles)
+    print(f"A_matrix shape: {A_matrix.shape}")
+
+    print("\nFirst 5 pole coordinates for C1:")
+    for i in range(5):
+        print(f"({poles_x1[i]:.2f}, {poles_y1[i]:.2f})")
+
+    print("\nFirst 5 pole coordinates for C2:")
+    for i in range(5):
+        print(f"({poles_x2[i]:.2f}, {poles_y2[i]:.2f})")
+
+    print("\nFirst 5 collocation points for C1:")
+    for i in range(5):
+        print(f"({coll_x1[i]:.2f}, {coll_y1[i]:.2f})")
+
+    print("\nFirst 5 collocation points for C2:")
+    for i in range(5):
+        print(f"({coll_x2[i]:.2f}, {coll_y2[i]:.2f})")
+
+    print("\nFirst 10 values of b_vector:")
+    for i in range(10):
+        print(f"{b_vector[i]:.2f}")
+
+    print("\nFirst 5x5 submatrix of A:")
+    # Ensure we only print up to 5x5, even if matrix is smaller
+    rows_to_print = min(5, A_matrix.shape[0])
+    cols_to_print = min(5, A_matrix.shape[1])
+    for i in range(rows_to_print):
+        row_str = " ".join([f"{A_matrix[i, j]:9.4f}" for j in range(cols_to_print)])
+        print(f"[{row_str} ]")