✨ V1

Co-authored-by: tomchauveau <[email protected]>

+tools/+utils/digCube.m

@@ -0,0 +1,32 @@
+function geometry = digCube(geometry)
+    count = 0; 
+    while count < geometry.nbr_cavities
+        radius = geometry.min_radius+(geometry.max_radius-geometry.min_radius)*rand();
+        x = (geometry.unit-2*radius)*rand(1)-geometry.unit/2+radius;
+        y = (geometry.unit-2*radius)*rand(1)-geometry.unit/2+radius;
+        z = (geometry.unit-2*radius)*rand(1)-geometry.unit/2+radius;
+        valid = true; 
+        for i=1:length(geometry.radii) 
+            if tools.utils.spheresIntersect([x,y,z],radius,geometry.centers(i,:),geometry.radii(i)) 
+                valid = false; 
+               break; 
+            end
+        end
+
+        if valid 
+            geometry.radii = [geometry.radii, radius]; 
+            geometry.centers = [geometry.centers; x, y, z]; 
+            cavity = multisphere(radius);
+            cavity = translate(cavity,[x y z]);
+            if strcmp(geometry.type, 'full')
+                geometry.structure = addCell(geometry.structure,cavity);
+            elseif strcmp(geometry.type, 'void')
+                geometry.structure = addVoid(geometry.structure,cavity);
+            end
+            count = count+1; 
+            geometry.porous_volume = geometry.porous_volume + 4/3*pi*radius^3;
+            geometry.porosity_fraction = (geometry.porous_volume/geometry.volume)*100;
+        end
+
+    end
+end

+tools/+utils/digEllipsoid.m

@@ -0,0 +1,45 @@
+function geometry = digEllipsoid(geometry)
+    count = 0;
+    while count < geometry.nbr_cavities
+        radius = geometry.min_radius+(geometry.max_radius-geometry.min_radius)*rand();       
+        valid = false;
+        while ~valid
+            x = 2*geometry.unit(1)*(rand(1)-0.5);
+            y = 2*geometry.unit(2)*(rand(1)-0.5);
+            z = 2*geometry.unit(3)*(rand(1)-0.5);
+            if isInsideEllipsoid(x, y, z, geometry.unit(1), geometry.unit(2), geometry.unit(3), radius)
+                valid = true;
+            end
+        end
+
+        for i = 1:length(geometry.radii)
+            if tools.utils.spheresIntersect([x, y, z], radius, geometry.centers(i,:), geometry.radii(i))
+                valid = false;
+                break;
+            end
+        end
+
+        if valid
+            geometry.radii = [geometry.radii, radius];
+            geometry.centers = [geometry.centers; x, y, z];
+
+            cavity = multisphere(radius);
+            cavity = translate(cavity, [x y z]);
+            if strcmp(geometry.type, 'full')
+                geometry.structure = addCell(geometry.structure,cavity);
+            elseif strcmp(geometry.type, 'void')
+                geometry.structure = addVoid(geometry.structure,cavity);
+            end
+            count = count+1;
+            geometry.porous_volume = geometry.porous_volume + 4/3*pi*radius^3;
+            geometry.porosity_fraction = (geometry.porous_volume/geometry.volume)*100;
+        end
+    end
+end
+
+function inside = isInsideEllipsoid(x, y, z, a, b, c, radius)
+    x_norm = abs(x) + radius;
+    y_norm = abs(y) + radius;
+    z_norm = abs(z) + radius;
+    inside = (x_norm / a)^2 + (y_norm / b)^2 + (z_norm / c)^2 <= 1;
+end

+tools/+utils/spheresIntersect.m

@@ -0,0 +1,4 @@
+function intersect = spheresIntersect(center1, radius1, center2, radius2)
+    distance = sqrt(sum((center1 - center2).^2));
+    intersect = distance < (radius1 + radius2);
+end

+tools/generateGeometry.m

@@ -0,0 +1,29 @@
+function geometry = generateGeometry(shape, unit, nbr_cavities, type, min_radius, max_radius)
+    geometry.unit = unit;
+    geometry.min_radius = min_radius;
+    geometry.max_radius = max_radius;
+    geometry.nbr_cavities = nbr_cavities;
+    if ismember(type, {'void', 'full'})
+        geometry.type = type;
+    else
+        error("Unknown geometry type: must be 'void' or 'full' cavities.");
+    end
+    geometry.centers = [0,0,0];
+    geometry.porous_volume = 0; 
+    geometry.porosity_fraction = 0;
+    geometry.radii = 0.1;
+
+    if strcmp(shape, 'cube')
+        geometry.volume = geometry.unit^3;
+        geometry.structure = multicuboid(geometry.unit,geometry.unit,geometry.unit,'Zoffset',-geometry.unit/2);
+        geometry = tools.utils.digCube(geometry);
+        geometry.exposed_faces=6;
+    elseif strcmp(shape, 'ellipsoid')
+        geometry.volume = 4/3*pi*unit(1)*unit(2)*unit(3);
+        geometry.structure = scale(multisphere(1),[unit(1),unit(2),unit(3)]);
+        geometry = tools.utils.digEllipsoid(geometry);
+        geometry.exposed_faces=1;
+    else
+        error("generateGeometry : Invalid shape argument")
+    end
+end

+tools/geometryView.m

@@ -0,0 +1,25 @@
+function geometryView(filename, thermalResults, geometry)
+    % Parameters for figure customization
+    set(groot,'defaultAxesTickLabelInterpreter','latex'); 
+    set(groot,'defaulttextinterpreter','latex');
+    set(groot,'defaultLegendInterpreter','latex');
+
+    tempFigure = figure('Visible', 'off');
+
+    pdeplot3D(thermalResults.Mesh,ColorMapData=thermalResults.Temperature(:,end));
+    title({sprintf("Solution plot at %0.2f days", thermalResults.SolutionTimes(end)/(24*3600)), " "}, 'FontSize', 15);
+    exportgraphics(tempFigure, filename, 'ContentType', 'vector', 'Resolution', 700); 
+
+    clf(tempFigure);
+
+    pdegplot(geometry.structure, "FaceAlpha", 0.2);
+    title({"Geometry representation", " "}, 'FontSize', 15);
+    exportgraphics(tempFigure, filename, 'ContentType', 'vector', 'Resolution', 700, 'Append', true); 
+
+    clf(tempFigure);
+
+    pdemesh(thermalResults.Mesh);
+    title({"Mesh representation", " "}, 'FontSize', 15);
+    exportgraphics(tempFigure, filename, 'ContentType', 'vector', 'Resolution', 700, 'Append', true); 
+
+end

+tools/getThermalBarycenter.m

@@ -0,0 +1,15 @@
+function Result = getThermalBarycenter(thermalResults, Temperature)
+    A=thermalResults.Temperature;
+    lastPointIndex=-1;
+    lastTimeIndex=-1;
+    for t=size(A,2):-1:1
+        pointIndices = find(A(:,t) >= Temperature);
+        if ~isempty(pointIndices)
+            lastPointIndex=pointIndices(1);
+            lastTimeIndex=t;
+            break;
+        end
+    end
+    Result.Bar = thermalResults.Mesh.Nodes(:,lastPointIndex);
+    Result.Time = thermalResults.SolutionTimes(1,lastTimeIndex);
+end

+tools/heatMapSlice.m

@@ -0,0 +1,69 @@
+function heatMapSlice(filename,thermalResults,view_type,value)
+
+% Parameters for figure customization
+set(groot,'defaultAxesTickLabelInterpreter','latex'); 
+set(groot,'defaulttextinterpreter','latex');
+set(groot,'defaultLegendInterpreter','latex');
+
+
+Mesh = thermalResults.Mesh;
+x = linspace(min(Mesh.Nodes(1,:)), max(Mesh.Nodes(1,:)), 200);
+y = linspace(min(Mesh.Nodes(2,:)), max(Mesh.Nodes(2,:)), 200);
+z = linspace(min(Mesh.Nodes(3,:)), max(Mesh.Nodes(3,:)), 200); 
+
+if view_type == 'x' 
+    [Y,Z] = meshgrid(y,z); 
+    X = value*ones(size(Y));
+elseif view_type == 'y'
+    [X,Z] = meshgrid(x,z); 
+    Y = value*ones(size(X));
+elseif view_type == 'z'
+    [X,Y] = meshgrid(x,y);
+    Z = value*ones(size(X));
+else 
+    error('Invalid view type. Use ''z'', ''x'', or ''y''.');
+end
+
+
+day = 24*3600;
+tIndex = [100, 300, 500, 700, 900, 1100];
+Temp = interpolateTemperature(thermalResults, X, Y, Z, tIndex);
+
+figure("Position",[10,10,1000,550], 'Visible', 'off'); 
+hold on; 
+for i = 1:numel(tIndex)
+    subplot(2, 3, i)
+
+    if view_type == 'x'
+        TempInterp = interp2(y, z, reshape(Temp(:, i), size(Y)), Y, Z, 'bicubic');
+        surf(Y, Z, TempInterp);
+        xlabel("$y$ $(m)$");
+        ylabel("$z$ $(m)$");
+    elseif view_type == 'y'
+        TempInterp = interp2(x, z, reshape(Temp(:, i), size(X)), X, Z, 'bicubic');
+        surf(X, Z, TempInterp);
+        xlabel("$x$ $(m)$");
+        ylabel("$z$ $(m)$");
+    elseif view_type == 'z'
+        TempInterp = interp2(x, y, reshape(Temp(:, i), size(X)), X, Y, 'bicubic');
+        surf(X, Y, TempInterp);
+        xlabel("$x$ $(m)$");
+        ylabel("$y$ $(m)$");
+    end
+
+    title(sprintf('t = %1.2f hours', thermalResults.SolutionTimes(tIndex(i))/3600));
+    view(0,90);
+    shading interp; 
+    colormap("jet");
+    grid off; 
+    clim([300,2000]);
+end
+c = colorbar('Position', [0.93, 0.11, 0.02, 0.815]);
+c.Label.Interpreter = 'latex';
+c.TickLabelInterpreter = 'latex';
+
+
+% Save figure in file
+set(gcf, 'PaperOrientation', 'landscape');
+print(gcf, '-dpdf', '-r500', '-bestfit', filename);
+end

+tools/heatMapSliceAnimation.m

@@ -0,0 +1,76 @@
+function heatMapSliceAnimation(filename,thermalResults,view_type,value)
+
+% Parameters for figure customization
+set(groot,'defaultAxesTickLabelInterpreter','latex'); 
+set(groot,'defaulttextinterpreter','latex');
+set(groot,'defaultLegendInterpreter','latex');
+
+
+Mesh = thermalResults.Mesh;
+x = linspace(min(Mesh.Nodes(1,:)), max(Mesh.Nodes(1,:)), 400);
+y = linspace(min(Mesh.Nodes(2,:)), max(Mesh.Nodes(2,:)), 400);
+z = linspace(min(Mesh.Nodes(3,:)), max(Mesh.Nodes(3,:)), 400); 
+
+if view_type == 'x' 
+    [Y,Z] = meshgrid(y,z); 
+    X = value*ones(size(Y));
+elseif view_type == 'y'
+    [X,Z] = meshgrid(x,z); 
+    Y = value*ones(size(X));
+elseif view_type == 'z'
+    [X,Y] = meshgrid(x,y);
+    Z = value*ones(size(X));
+else 
+    error('Invalid view type. Use ''z'', ''x'', or ''y''.');
+end
+
+
+day = 24*3600;
+timestep = 5; % to accelerate rendering
+Temp = interpolateTemperature(thermalResults, X, Y, Z, [1:size(thermalResults.SolutionTimes,2)]);
+
+figure('Visible', 'off'); 
+hold on; 
+video = VideoWriter(filename);
+video.FrameRate = 10; 
+open(video);
+for i = 1:timestep:size(thermalResults.SolutionTimes,2) 
+    if view_type == 'x'
+        TempInterp = interp2(y, z, reshape(Temp(:, i), size(Y)), Y, Z, 'bicubic');
+        p = surf(Y, Z, TempInterp);
+        xlabel("$y$ $(m)$");
+        ylabel("$z$ $(m)$");
+    elseif view_type == 'y'
+        TempInterp = interp2(x, z, reshape(Temp(:, i), size(X)), X, Z, 'bicubic');
+        p = surf(X, Z, TempInterp);
+        xlabel("$x$ $(m)$");
+        ylabel("$z$ $(m)$");
+    elseif view_type == 'z'
+        TempInterp = interp2(x, y, reshape(Temp(:, i), size(X)), X, Y, 'bicubic');
+        p = surf(X, Y, TempInterp);
+        xlabel("$x$ $(m)$");
+        ylabel("$y$ $(m)$");
+    end
+
+    title(sprintf('t = %1.2f hours\n', thermalResults.SolutionTimes(i)/3600), 'FontSize', 16);
+    view(0,90);
+    shading interp; 
+    colormap("jet");
+    grid off; 
+    clim([300,2000]);
+    c = colorbar; 
+    c.Label.Interpreter = 'latex';
+    c.TickLabelInterpreter = 'latex';
+    ylabel(c, 'T (K)', 'FontSize', 10, 'Interpreter', 'latex');
+    c.Label.Position(1) = c.Label.Position(1) + 1;
+
+    pause(0.01);
+    frame = getframe(gcf);
+    writeVideo(video, frame);
+    if i < size(thermalResults.SolutionTimes,2)
+        delete(p);
+    end
+
+end
+close(video);
+end

+tools/isosurface.m

@@ -0,0 +1,42 @@
+function isosurface(filepath, thermalResults, temperature)
+    %% Loading Data
+
+    Mesh = thermalResults.Mesh; 
+
+    %% Animation
+    figure('Visible','off');
+
+    % Export to video file
+    video = VideoWriter(filepath);
+    video.FrameRate = 10; 
+    open(video);
+
+    hold on;
+    axis off
+    pdeplot3D(Mesh.Nodes, Mesh.Elements, 'FaceAlpha', 0.1, 'EdgeColor', 'none');
+    for t = 1:size(thermalResults.SolutionTimes,2)
+        temp = thermalResults.Temperature(:, t);
+        F = scatteredInterpolant(Mesh.Nodes(1,:)', Mesh.Nodes(2,:)', Mesh.Nodes(3,:)', temp, 'linear', 'none');
+        [x, y, z] = meshgrid(linspace(min(Mesh.Nodes(1,:)), max(Mesh.Nodes(1,:)), 50), ...
+            linspace(min(Mesh.Nodes(2,:)), max(Mesh.Nodes(2,:)), 50), ...
+                linspace(min(Mesh.Nodes(3,:)), max(Mesh.Nodes(3,:)), 50));
+
+                v = F(x, y, z);
+                axis equal 
+        %axis vis3d  % rotate the 3D animation
+        p = patch(isosurface(x, y, z, v, temperature));
+        p.FaceColor = 'blue';
+        p.EdgeColor = 'none'; 
+        p.FaceAlpha = 0.2; % opacity
+        view(3)
+        colormap copper
+        camlight; lighting phong;
+        pause(0.1);
+        frame = getframe(gcf);
+        writeVideo(video, frame);
+        if t < size(thermalResults.SolutionTimes,2)
+            delete(p);
+        end
+    end
+    close(video);
+    hold off;

+tools/modifyGeometry.m

@@ -0,0 +1,25 @@
+function geometry = modifyGeometry(previous_geometry, shape, unit, nbr_cavities, type, min_radius, max_radius)
+    geometry.structure = previous_geometry;
+    geometry.unit = unit;
+    geometry.min_radius = min_radius;
+    geometry.max_radius = max_radius;
+    geometry.nbr_cavities = nbr_cavities;
+    if ismember(type, {'void', 'full'})
+        geometry.type = type;
+    else
+        error("Unknown geometry type: must be 'void' or 'full' cavities.");
+    end
+    geometry.centers = [0,0,0];
+    geometry.porous_volume = 0; 
+    geometry.porosity_fraction = 0;
+    geometry.radii = 0.1;
+    if strcmp(shape, 'cube')
+        geometry.volume = geometry.unit^3;
+        geometry = tools.utils.digCube(geometry);
+    elseif strcmp(shape, 'ellipsoid')
+        geometry.volume = 4/3*pi*unit(1)*unit(2)*unit(3);
+        geometry = tools.utils.digEllipsoid(geometry);
+    else
+        error("modifyGeometry : Invalid shape argument")
+    end
+end

.gitignore

@@ -0,0 +1,17 @@
+.DS_Store
+__pycache__
+simulations
+*.mat
+*.stl
+*.avi
+*.log
+*.fls
+*.out
+*.gz
+*.toc
+*.fdb_latexmk
+*.aux
+script.m
+script.py
+simulations.ipynb
+doc/latex-UG/main.pdf