% MATLAB Startup preparation
close all; clc; clear; format compact;
% Define Mesh Parameters for X, Y, and Z axes
expansionFactor = 1.1;
numCellsLinearX = 12; numCellsLinearY = 12; numCellsLinearZ = 8;
startLinearX = 0; endLinearX = 4.6;
startLinearY = 0; endLinearY = 4.6;
startLinearZ = 0; endLinearZ = 3;
startNonlinearBeforeX = -9.2; endNonlinearAfterX = 18.4;
startNonlinearBeforeY = -9.2; endNonlinearAfterY = 13.8;
startNonlinearBeforeZ = 0; endNonlinearAfterZ = 9;
% Generate points for X, Y, and Z axes
xPoints = generatePoints(startLinearX, endLinearX, startNonlinearBeforeX, endNonlinearAfterX, numCellsLinearX, expansionFactor);
yPoints = generatePoints(startLinearY, endLinearY, startNonlinearBeforeY, endNonlinearAfterY, numCellsLinearY, expansionFactor);
zPoints = generatePoints(startLinearZ, endLinearZ, startNonlinearBeforeZ, endNonlinearAfterZ, numCellsLinearZ, expansionFactor);
% Create mesh grid
[X, Y, Z] = meshgrid(xPoints, yPoints, zPoints);
% Calculate the volume of each cell (assuming rectangular prisms)
cellVolumes = diff(xPoints) .* diff(yPoints)' .* permute(diff(zPoints), [3, 1, 2]);
cellVolumes = cellVolumes(:);
% Identify indices for linear and nonlinear regions
linearMask = X(1:end-1, 1:end-1, 1:end-1) >= startLinearX & X(1:end-1, 1:end-1, 1:end-1) <= endLinearX & ...
Y(1:end-1, 1:end-1, 1:end-1) >= startLinearY & Y(1:end-1, 1:end-1, 1:end-1) <= endLinearY & ...
Z(1:end-1, 1:end-1, 1:end-1) >= startLinearZ & Z(1:end-1, 1:end-1, 1:end-1) <= endLinearZ;
linearVolumes = cellVolumes(linearMask);
% Calculate the characteristic length (h) for linear and nonlinear regions
h_linear = mean(nthroot(linearVolumes, 3));
h_nonlinear = mean(nthroot(cellVolumes(~linearMask), 3));
h_total = mean(nthroot(cellVolumes, 3));
% Calculate the number of cells in each region
numCellsLinear = numel(linearVolumes);
numCellsNonlinear = numel(cellVolumes) - numCellsLinear;
numCellsTotal = numel(cellVolumes);
% Prepare the table data
cellCounts = [numCellsLinear; numCellsNonlinear; numCellsTotal];
regionNames = {'Linear Region', 'Nonlinear Region', 'Total'};
hValues = [h_linear; h_nonlinear; h_total]; % Characteristic lengths for each region
% Create and display the table
cellCountTable = table(regionNames, cellCounts, hValues, 'VariableNames', {'Region', 'NumberOfCells', 'CharacteristicLength'});
disp(cellCountTable);
% Visualization
visualizeMesh(X, Y, Z, linearMask);
% Function to generate linear and nonlinear points
function points = generatePoints(startLinear, endLinear, startNonlinear, endNonlinear, numCellsLinear, expansionFactor)
cellSizeLinear = (endLinear - startLinear) / numCellsLinear;
pointsLinear = linspace(startLinear, endLinear, numCellsLinear + 1);
pointsNonlinearBefore = generateNonlinearPoints(startLinear, -cellSizeLinear, expansionFactor, true, startNonlinear);
pointsNonlinearAfter = generateNonlinearPoints(endLinear, cellSizeLinear, expansionFactor, false, endNonlinear);
points = sort(unique([pointsNonlinearBefore, pointsLinear, pointsNonlinearAfter]));
end
% Visualization function
function visualizeMesh(X, Y, Z, linearMask)
figure;
hold on;
plot3(X(linearMask), Y(linearMask), Z(linearMask), 'r.', 'MarkerSize', 10);
plot3(X(~linearMask), Y(~linearMask), Z(~linearMask), 'b.', 'MarkerSize', 10);
xlabel('X-axis'); ylabel('Y-axis'); zlabel('Z-axis');
title('3D Mesh with Linear and Nonlinear Regions');
axis equal; grid on; view(3);
hold off;
end