Multi-body

.gitignore

@@ -34,4 +34,13 @@ octave-workspace
 *.log
 
 # Virtual environment
-venv/
+venv/
+
+# Images
+*.png
+*.jpg
+*.jpeg
+*.gif
+*.bmp
+*.svg
+*.eps

config.m

@@ -10,9 +10,10 @@
     %   cfg = config('ellipse');           % Ellipse geometry
     %   cfg = config('rectangle');         % Rectangle geometry
     %   cfg = config('airfoil');           % NACA airfoil geometry
+    %   cfg = config('multi');             % Multiple obstacles geometry
 
     if nargin < 1
-        geometry_type = 'rectangle';  % Default geometry
+        geometry_type = 'multi';  % Default geometry
     end
 
     %% Domain Configuration
@@ -30,10 +31,21 @@
     config.mesh.refine_b2 = 0.08;             % Mesh refinement parameter B2 (wake region)
     config.mesh.edge_multiplier = 3;          % Multiplier for edge generation from triangles
 
+    %% Simulation Parameters (set defaults first, then override per geometry)
+    config.simulation.reynolds_number = 100;
+    config.simulation.viscosity = 1 / config.simulation.reynolds_number;
+    config.simulation.time_step = 1e-2;  % Default time step (may be overridden by geometry)
+    config.simulation.num_time_steps = 2000;
+    config.simulation.num_time_steps_ci = 20;
+    config.simulation.random_seed = 42;  % Required for DistMesh reproducibility (uses rand() for rejection method)
+    config.simulation.show_progress = true;  % Display time step progress (disabled in CI)
+
     %% Geometry Parameters - Set based on geometry type
-    config.geometry.type = lower(geometry_type);
+    % Normalize geometry type to lowercase once
+    geometry_type_lower = lower(geometry_type);
+    config.geometry.type = geometry_type_lower;
 
-    switch lower(geometry_type)
+    switch geometry_type_lower
         case 'cylinder'
             config.geometry.obstacle_radius = 0.5;    % Cylinder radius
 
@@ -48,18 +60,31 @@
             config.geometry.rect_y_center = 0.0;      % Rectangle center Y-coordinate
             % Rectangle needs finer mesh for convergence
             config.mesh.dist = 0.05;
+
         case 'airfoil'
             % NACA 4-digit series parameters
             config.geometry.naca_digits = [0, 0, 1, 9];       % NACA airfoil
             config.geometry.chord_length = 2.3;               % Chord length
             config.geometry.angle_of_attack = -10;            % Angle of attack in degrees
             config.geometry.airfoil_x_center = -0.5;          % Airfoil center X-coordinate (leading edge, shifted left)
             config.geometry.airfoil_y_center = 0.0;           % Airfoil center Y-coordinate
+
+            config.mesh.dist = 0.075;                          % Need a somewhat finer mesh distribution
+            config.simulation.time_step = 5e-3;               % Airfoil needs smaller time step for stability
+
+        case 'multi'
+            % Multiple obstacles configuration
+            % Default: two cylinders side-by-side vertically
+            config.geometry.obstacles = [ ...
+                                         struct('type', 'cylinder', 'center', [0, 2.5], 'params', struct('radius', 0.5)), ...
+                                         struct('type', 'cylinder', 'center', [0, -2.5], 'params', struct('radius', 0.5)) ...
+                                        ];
+
             config.mesh.dist = 0.05;                          % Need a somewhat finer mesh distribution
-            config.mesh.refine_a1 = 0.02;                     % Mesh refinement parameter A1 (near obstacle)
-            config.mesh.refine_b1 = 0.05;                     % Mesh refinement parameter B1 (near obstacle)
+            config.simulation.time_step = 5e-3;
+
         otherwise
-            error('Unknown geometry type: %s. Supported types: cylinder, ellipse, rectangle, airfoil', geometry_type);
+            error('Unknown geometry type: %s. Supported types: cylinder, ellipse, rectangle, airfoil, multi', geometry_type);
     end
 
     %% RBF-FD Algorithm Parameters
@@ -92,21 +117,7 @@
     config.rbf.order_near_boundary = 17;
     config.rbf.poly_degree_near_boundary = 2;
 
-    %% Simulation Parameters
-    config.simulation.reynolds_number = 100;
-    config.simulation.viscosity = 1 / config.simulation.reynolds_number;
-    config.simulation.time_step = 1e-2;
-    config.simulation.num_time_steps = 1000;
-    config.simulation.num_time_steps_ci = 20;
-    config.simulation.random_seed = 42;  % Required for DistMesh reproducibility (uses rand() for rejection method)
-    config.simulation.show_progress = true;  % Display time step progress (disabled in CI)
-
-    % Geometry-specific time step adjustments for stability
-    if strcmp(lower(geometry_type), 'airfoil')
-        config.simulation.time_step = 5e-3;  % Smaller time step for airfoil stability
-    end
-
-    %% Distance Thresholds for Special Treatment
+    %% Boundary Distance Thresholds for Special Treatment
     config.distances.x_min = 1;
     config.distances.x_max = 1;
     config.distances.y_min = 0.5;

simulate.m

@@ -80,6 +80,27 @@
     thickness_percent = naca_digits(3) * 10 + naca_digits(4); % Thickness as percentage
     max_thickness = (thickness_percent / 100) * chord_length;  % Maximum thickness in units
     obs_char_length = max_thickness / 2; % Use half of maximum thickness
+elseif isfield(G, 'obstacles') && isfield(G, 'num_obstacles')
+    % Multi-obstacle geometry: use minimum characteristic length of all obstacles
+    obs_char_length = inf;
+    for i = 1:G.num_obstacles
+        obs = G.obstacles(i);
+        switch lower(obs.type)
+            case 'cylinder'
+                char_len = obs.params.radius;
+            case 'ellipse'
+                char_len = min(obs.params.a, obs.params.b);
+            case 'rectangle'
+                char_len = min(obs.params.width, obs.params.height) / 2;
+            case 'airfoil'
+                thickness_percent = obs.params.naca_digits(3) * 10 + obs.params.naca_digits(4);
+                max_thickness = (thickness_percent / 100) * obs.params.chord_length;
+                char_len = max_thickness / 2;
+            otherwise
+                char_len = 0.5; % Default fallback
+        end
+        obs_char_length = min(obs_char_length, char_len);
+    end
 else
     error('Geometry structure missing expected parameters');
 end
@@ -88,17 +109,91 @@
 boundary_s = [boundary_y_s; boundary_out_s; boundary_in_s; boundary_obs_s];
 boundary = [boundary_y; boundary_out; boundary_in; boundary_obs];
 
-% Initial mesh visualization
-if doPlot
-    figure;
-    scatter(xy(:, 1), xy(:, 2), 'k.');
+% Initial mesh visualization (only in non-CI mode)
+if doPlot && ~isCI && ~isTest
+    fprintf('Creating mesh visualization...\n');
+    fig = figure('Name', 'Multi-Obstacle Mesh', 'Position', [100, 100, 1000, 600]);
+
+    % Plot interior nodes (sample for performance if needed)
+    if size(xy, 1) > 10000
+        % Sample for very large meshes
+        sample_rate = max(1, floor(size(xy, 1) / 8000));
+        scatter(xy(1:sample_rate:end, 1), xy(1:sample_rate:end, 2), 4, 'k.', 'DisplayName', ...
+                sprintf('Interior nodes (sampled %d/%d)', length(1:sample_rate:size(xy, 1)), size(xy, 1)));
+    else
+        % Show all interior nodes for reasonable mesh sizes
+        scatter(xy(:, 1), xy(:, 2), 3, 'k.', 'DisplayName', sprintf('Interior nodes (%d)', size(xy, 1)));
+    end
     hold on;
-    axis square;
-    scatter(boundary_in(:, 1), boundary_in(:, 2), 'b+');
-    scatter(boundary_y(:, 1), boundary_y(:, 2), 'r+');
-    scatter(boundary_out(:, 1), boundary_out(:, 2), 'b+');
-    scatter(boundary_obs(:, 1), boundary_obs(:, 2), 'm+');
+
+    % Plot pressure grid interior nodes (lighter color)
+    if size(xy_s, 1) > 5000
+        % Sample pressure nodes for very large meshes
+        p_sample_rate = max(1, floor(size(xy_s, 1) / 4000));
+        scatter(xy_s(1:p_sample_rate:end, 1), xy_s(1:p_sample_rate:end, 2), 2, [0.7 0.7 0.7], '.', ...
+                'DisplayName', sprintf('Pressure nodes (sampled %d/%d)', length(1:p_sample_rate:size(xy_s, 1)), size(xy_s, 1)));
+    else
+        scatter(xy_s(:, 1), xy_s(:, 2), 2, [0.7 0.7 0.7], '.', ...
+                'DisplayName', sprintf('Pressure nodes (%d)', size(xy_s, 1)));
+    end
+
+    % Plot boundary nodes with different colors and markers
+    scatter(boundary_in(:, 1), boundary_in(:, 2), 25, 'b', 's', 'filled', 'DisplayName', sprintf('Inlet (%d)', length(boundary_in)));
+    scatter(boundary_y(:, 1), boundary_y(:, 2), 25, 'r', '^', 'filled', 'DisplayName', sprintf('Walls (%d)', length(boundary_y)));
+    scatter(boundary_out(:, 1), boundary_out(:, 2), 25, 'g', 'v', 'filled', 'DisplayName', sprintf('Outlet (%d)', length(boundary_out)));
+
+    % Color-code obstacle boundaries if multi-obstacle
+    if isfield(G, 'boundary_obs_ids') && isfield(G, 'num_obstacles')
+        colors = {[1 0 1], [0 1 1], [1 1 0], [0.8 0.4 0], [0.6 0.2 0.8], [0.2 0.8 0.2], [0.8 0.2 0.2]};
+        for obs_id = 1:G.num_obstacles
+            obs_mask = G.boundary_obs_ids == obs_id;
+            if any(obs_mask)
+                color_idx = mod(obs_id - 1, length(colors)) + 1;
+                color = colors{color_idx};
+                scatter(G.boundary_obs_s(obs_mask, 1), G.boundary_obs_s(obs_mask, 2), ...
+                        35, color, 'o', 'filled', 'LineWidth', 1, 'MarkerEdgeColor', 'k', ...
+                        'DisplayName', sprintf('Obstacle %d (%d pts)', obs_id, sum(obs_mask)));
+            end
+        end
+    else
+        scatter(boundary_obs(:, 1), boundary_obs(:, 2), 35, [1 0 1], 'o', 'filled', ...
+                'LineWidth', 1, 'MarkerEdgeColor', 'k', ...
+                'DisplayName', sprintf('Obstacle (%d pts)', length(boundary_obs)));
+    end
+
     axis equal;
+    grid on;
+
+    % Calculate total node counts
+    total_v_nodes = size(xy, 1) + length(boundary_in) + length(boundary_y) + length(boundary_out) + length(boundary_obs);
+    total_p_nodes = size(xy_s, 1) + length(G.boundary_y_s) + length(G.boundary_out_s) + length(G.boundary_in_s) + length(G.boundary_obs_s);
+    if isfield(G, 'num_obstacles')
+        num_obstacles = G.num_obstacles;
+    else
+        num_obstacles = 1; % Single obstacle case
+    end
+
+    title(sprintf('Multi-Obstacle Mesh: %d obstacles, V-grid: %d nodes, P-grid: %d nodes', ...
+                  num_obstacles, total_v_nodes, total_p_nodes));
+    xlabel('x');
+    ylabel('y');
+    legend('Location', 'best', 'FontSize', 8);
+
+    % Save figure and display
+    try
+        saveas(fig, 'multi_obstacle_mesh.png');
+        fprintf('Mesh visualization saved as multi_obstacle_mesh.png\n');
+    catch
+        % Ignore save errors
+    end
+
+    % Force figure to front and pause briefly to ensure visibility
+    figure(fig);
+    drawnow;
+    fprintf('Mesh visualization displayed. Continuing with simulation...\n');
+    pause(1); % Brief pause to let user see the mesh
+elseif isCI || isTest
+    fprintf('Skipping mesh visualization (CI/test mode)\n');
 end
 
 % Combine interior and boundary nodes to create complete grids

src/build_geometry.m

@@ -31,8 +31,11 @@
                     cfg.geometry.chord_length, cfg.geometry.angle_of_attack, ...
                     cfg.geometry.airfoil_x_center, cfg.geometry.airfoil_y_center);
             G = make_airfoil_geometry(cfg);
+        case 'multi'
+            fprintf('Using multi-obstacle geometry (%d obstacles)...\n', length(cfg.geometry.obstacles));
+            G = make_multi_geometry(cfg);
         otherwise
-            error('Unknown geometry type: %s. Supported types: cylinder, ellipse, rectangle, airfoil', cfg.geometry.type);
+            error('Unknown geometry type: %s. Supported types: cylinder, ellipse, rectangle, airfoil, multi', cfg.geometry.type);
     end
 
 end