How can I use the isocaps function or any other method to close the surfaces of a cylindrical gyroid?

26 次查看（过去 30 天）
CHEICK IBRAHIM
CHEICK IBRAHIM 2024-11-8，14:28
回答： CHEICK IBRAHIM about 9 hours 前
Hello,
As part of my thesis, I've been asked to create a Matlab program to generate a cylindrical gyroid (this is essential in view of the dynamic test equipment to be used later). To do this, I was asked to start with a classic prismatic gyroid and then apply a cylindricity condition to generate a cylindrical model.
Problem: when I associated isosurface and isocaps to generate the basic prismatic gyroid to apply the cylindricity condition (I'll spare you the details.), I was getting errors on my radial surfaces.
So we proceeded without the isocaps function (code below) to create the cylindrical gyroid from the prismatic model. My current problem is to be able to close the surfaces as the isocaps function would.
%% cleanup section
close all
clear
clc
% Define envelope function
dimensionX = 10; % structure dimension in the x-direction
dimensionY = 10; % structure dimension in the y-direction
dimensionZ = 10; % structure dimension the z-direction
resolution = 0.5; % shape resolution (isotropic)
cellLength = 10; % the length of one gyroid cell
thickness = 1; % thickness of the gyroid surface
FV=createGyroid(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness);
f = @(x, y, z) z - 3;
createIntersectedVolume(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f,FV);
f = @(x, y, z) x.*x + y.*y -25;
createIntersectedVolume(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f,FV);
%%%%% Different functions used in these codes %%%%%%%
%% =========================================================================
function createIntersectedVolume(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f,FV)
%figure;
%plotStructure(FV);
%plotf(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f)
%% Search points outside envelope
xv = FV.vertices(:, 1);
yv = FV.vertices(:, 2);
zv = FV.vertices(:, 3);
fv=f(xv,yv,zv);
removePoints = find(fv > 0);
keepPoints = find(fv <= 0);
insideV.faces = FV.faces;
insideV.vertices = FV.vertices;
outsideV.faces = FV.faces;
outsideV.vertices = FV.vertices;
frontierV.faces = FV.faces;
frontierV.vertices= FV.vertices;
% Search points from outfind with a connection with points in infind
j=1;
%Remove faces whose points are outside
for i = size(FV.faces,1):-1:1
inpt1=find(FV.faces(i,1)==keepPoints);
inpt2=find(FV.faces(i,2)==keepPoints);
inpt3=find(FV.faces(i,3)==keepPoints);
%All points are outside
if (isempty(inpt1) && isempty(inpt2) && isempty(inpt3))
FV.faces(i,:)=[];
insideV.faces(i,:)=[];
frontierV.faces(i,:)=[];
%All points are inside
elseif ~(isempty(inpt1) || isempty(inpt2) || isempty(inpt3))
outsideV.faces(i,:)=[];
frontierV.faces(i,:)=[];
else
insideV.faces(i,:)=[];
outsideV.faces(i,:)=[];
end
end
for i = size(FV.faces,1):-1:1
ptin=[];
ptout=[];
for j=1:3
inptj=find(FV.faces(i,j)==keepPoints);
if ~isempty(inptj)
ptin=[ptin FV.faces(i,j)];
else
ptout=[ptout FV.faces(i,j)];
end
end
%% 2 point are inside, 1 point outside
tolerance = 10e-2; % Tolerance for comparing vertices
if size(ptout,2)==1
if i==1108
disp('');
end
%% Inspection of cases where the vertex already exists or does not exist after modification of his coordinates
[xNew, yNew, zNew] = fintersect(FV.vertices(ptin(1),1),FV.vertices(ptin(1),2),FV.vertices(ptin(1),3),...
FV.vertices(ptout,1), FV.vertices(ptout,2), FV.vertices(ptout,3),f);
newCoords=[xNew,yNew,zNew];
% Search for the Index of the vertex with the same coordinates after modification
existingIndex = findExistingVertex(FV.vertices, newCoords, tolerance);
if ~isempty(existingIndex)
% The modified vertex is similar to an existing vertex
% Update faces to use existing vertex index
FV.faces(FV.faces == ptout) = existingIndex;
ptId=existingIndex;
else
% Update vertex coordinates directly
FV.vertices(ptout, [1,2,3]) = newCoords;
ptId=ptout;
end
%% Add a point and a face because intersectin create 2 points
% Additional vertex created, regardless of whether a similar vertex exists
ind=size(FV.vertices,1)+1;
[FV.vertices(ind,1), FV.vertices(ind,2), FV.vertices(ind,3)] = ...
fintersect(FV.vertices(ptin(2),1),FV.vertices(ptin(2),2),FV.vertices(ptin(2),3),...
FV.vertices(ptout,1), FV.vertices(ptout,2), FV.vertices(ptout,3),f);
% Inspection of cases where the vertex already exists or does not exist
% Retrieve additional vertex coordinates
newVertex=[FV.vertices(ind,1), FV.vertices(ind,2), FV.vertices(ind,3)];
%Search for the Index of the vertex with the same coordinates as the one created
existingIndex = findExistingVertex(FV.vertices, newVertex, tolerance);
if isempty(existingIndex)
% Keep the new vertex and add a new face
FV.faces(end+1,1:3)=[ptin(2) ptId ind];
else
% Uses existing index if vertex already present
% FV.vertices(ind,:)=[]; % Delete the new vertex
ind=existingIndex;
FV.faces(end+1,1:3)=[ptin(2) ptId ind];
end
%% One point is inside, 2 points outside
elseif size(ptout,2)==2
for j=1:2
[xNew, yNew, zNew] = fintersect(FV.vertices(ptin,1),FV.vertices(ptin,2),FV.vertices(ptin,3),...
FV.vertices(ptout(j),1), FV.vertices(ptout(j),2), FV.vertices(ptout(j),3),f);
newCoords=[xNew,yNew,zNew];
% Search for the Index of the vertex with the same coordinates after modification
existingIndex = findExistingVertex(FV.vertices, newCoords, tolerance);
if ~isempty(existingIndex)
% The modified vertex is similar to an existing vertex
% Update faces to use existing vertex index
FV.faces(FV.faces == ptout(j)) = existingIndex;
else
% Update vertex coordinates directly
FV.vertices(ptout(j), :) = newCoords;
end
end
elseif ~isempty(ptout)
erreur('cas non prévu');
end
end
% Verification step
if any(FV.faces(:) > size(FV.vertices, 1))
error('Some indices in FV.faces exceed the size of FV.vertices');
end
%Remove unused vertices (i.e vertices not referenced by any face)
Vertices_to_keep=false(size(FV.vertices,1),1);
for ii=1:size(FV.faces,1)
Vertices_to_keep(FV.faces(ii,:))=true;
end
%Update the vertices
FV.vertices=FV.vertices(Vertices_to_keep,:);
%Remap the vertices in FV.faces
mapOld_toNew=zeros(size(Vertices_to_keep,1));
mapOld_toNew(Vertices_to_keep)=1:sum(Vertices_to_keep);
for ii=1:size(FV.faces,1)
FV.faces(ii,:)=mapOld_toNew(FV.faces(ii,:));
end
figure;
plotStructure(FV);
plotf(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f)
end
%% =========================================================================
function FV=createGyroid(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness)
% create evenly spaced points in each direction
% ---------------------------------------------
x = linspace(-dimensionX/2, dimensionX/2, dimensionX/resolution);
y = linspace(-dimensionY/2, dimensionY/2, dimensionY/resolution);
z = linspace(-dimensionZ/2, dimensionZ/2, dimensionZ/resolution);
% create a 3D grid from x, y and z
% --------------------------------
[X, Y, Z] = meshgrid(x, y, z);
% compute each point of two gyroids
% (upper/lower shells), and their intersection
% --------------------------------------------
F1 = cos(2*pi*X/cellLength).*sin(2*pi*Y/cellLength) + ...
cos(2*pi*Y/cellLength).*sin(2*pi*Z/cellLength) + ...
cos(2*pi*Z/cellLength).*sin(2*pi*X/cellLength) + ...
thickness/2;
F2 = cos(2*pi*X/cellLength).*sin(2*pi*Y/cellLength) + ...
cos(2*pi*Y/cellLength).*sin(2*pi*Z/cellLength) + ...
cos(2*pi*Z/cellLength).*sin(2*pi*X/cellLength) - ...
thickness/2;
F3 = F1.*F2;
% split the grid based on the expected
% walls thickness and close the domain
% ------------------------------------
FV = isosurface(X, Y, Z, F3, 0);
end
%% =========================================================================
function plotf(dimensionX,dimensionY,dimensionZ,resolution,cellLength,thickness,f)
x = linspace(-dimensionX/2, dimensionX/2, dimensionX/resolution);
y = linspace(-dimensionY/2, dimensionY/2, dimensionY/resolution);
z = linspace(-dimensionZ/2, dimensionZ/2, dimensionZ/resolution);
[X, Y, Z] = meshgrid(x, y, z);
F = f(X, Y, Z);
hold on;
p1 = patch(isosurface(X, Y, Z, F, 0));
p1.FaceColor = 'blue';
p1.EdgeColor = 'none';
p1.FaceAlpha = 0.2; % Transparency
hold off;
end
%% =========================================================================
function plotStructure(FV)
% create, plot and setup the plot
% of the patch created from the FV struct
% ---------------------------------------
% plot
pFV = patch(FV);
face_colors = rand(size(FV.faces, 1), 3); % Couleurs aléatoires en RGB% set style
%set(pFV, 'FaceColor', 'red', 'EdgeColor', 'none');
set(pFV, 'FaceVertexCData', face_colors, 'FaceColor', 'flat', 'EdgeColor', 'none');
hold off;
% set view and view point
daspect([1 1 1])
view([120, 15])
% set lighting
camlight
lighting flat
xlabel('X')
ylabel('Y')
zlabel('Z')
end
%% =========================================================================
function [x, y, z] = fintersect(x1, y1, z1, x2, y2, z2, f)
%REmarque : plutot que des fonctions anonymes, on pourrait appeler des
%fonctions.
%Par exemple, xt = @(t) x1 + t * (x2 - x1); on pourrait appeler une
%fonction avec les paramètres d'entrée x1,x2, t et en sortie xt la valeur
% de x entre x1 et x2 pour une valeur de t fixé
% Paramétrization x(t),t y(t), z(t) with t parameter in [0,1]
xt = @(t) x1 + t * (x2 - x1);
yt = @(t) y1 + t * (y2 - y1);
zt = @(t) z1 + t * (z2 - z1);
% surface with parametrization f(t)
tsurface = @(t) f(xt(t), yt(t), zt(t));
% Solve f(t)=0, intial value t=0.5
tsol = fzero(tsurface, 0.5);
x = xt(tsol);
y = yt(tsol);
z = zt(tsol);
end
function index = findExistingVertex(vertices, newVertex, tolerance)
% Cette fonction vérifie si un sommet similaire à `newVertex` existe déjà
% dans `vertices` avec une tolérance donnée. Elle retourne l'indice
% du premier sommet trouvé, sinon elle retourne [].
% Calcule la distance euclidienne entre `newVertex` et tous les sommets existants
diffs = vecnorm(vertices - newVertex, 2, 2);
% Trouve le premier sommet dont la distance est inférieure à la tolérance
index = find(diffs < tolerance, 1);
end
Sorry for the length of the code, but is there anyone who can help me solve my problem because I've been stuck on it for two weeks.
Thanks!!!!

请先登录，再回答此问题。

回答（2 个）

TED MOSBY
TED MOSBY about 12 hours 前
编辑：TED MOSBY about 12 hours 前
Hi @CHEICK IBRAHIM,
I ran your code that produces a gyroid using ‘isosurface’, the output looks like this:
You can use ‘isocaps’ for this problem as you mentioned but as you didn’t attach the code and specific issues you are facing using ‘isocaps’ I am attaching an example workflow on how to use the function such that your gyroid looks like this:
If you want to achieve something like this image, here is how I modified your code using ‘isocaps’ to achieve that output:
Refer to the code file 'code.m' in the attachments.
The new function ‘createCylindricalGyroid’ uses isocaps, which automatically adds caps at the open ends of the surface. The gyroid surface is plotted using patch(FV) and visualized with isocaps(X, Y, Z, F3, 0), ensuring that the open ends of the gyroid are automatically closed by adding caps.
The ‘createIntersectedVolume’ function in the code involves a cylindrical trim. It checks whether the points are inside or outside the cylindrical boundary defined by the function f(x, y, z) = x^2 + y^2 - radius^2, and then selectively retains or removes the faces and vertices that are inside the cylinder. It also handles the modification of faces for intersections.
Please refer to this MATLAB answer if you want to cap your structure like here: www.mathworks.com/support/search.html/answers/2032484-how-do-i-add-isocaps-to-the-other-voids-around-a-gyroid.html?fq%5B%5D=asset_type_name:answer&fq%5B%5D=category:matlab/volume-visualization&page=1
Hope this helps!
CHEICK IBRAHIM
CHEICK IBRAHIM about 9 hours 前
Hi @TED MOSBY
I thank you for the time you have dedicated to my problem.
Indeed, I deliberately omitted mentioning the issue I’m facing with isocaps because the result obtained was not the one I expected. To summarize, I create a gyroid by applying isocaps and then project all points outside the cylindrical domain, as shown in the attached screenshot
As can be seen, keeping all the faces by projecting the vertices does not produce the correct result. This is why I considered removing the faces that do not belong to the cylindrical domain. Here is the result I obtained (in the code that allowed me to achieve this, I also use isocaps before removing the faces):
As shown in this model, the removed faces are indeed missing. However, as I used isocaps to start with, the top and bottom faces are correct because the Z direction is not affected by the cylindrical condition imposed.
This is why, in the code I provided, I decided not to use isocaps at the beginning but instead see if I could apply it after obtaining the cylinder, which led to my question on the forum as I had exhausted all my ideas.
Regarding your suggestion, I reviewed your code and noticed that your function @createIntersectedVolume, which you call in createCylindricalGyroid to create the cylindrical gyroid, does not return anything. It therefore does not contribute to the creation of your cylindrical gyroid. Even when I change the radius, I still end up with a prism. Moreover, you directly apply isocaps on the base 3D matrices X, Y, Z, and F3.
This is exactly my issue: how to return the 3D matrices X, Y, Z, and F3 correctly after the function @createIntersectedVolume, in order to close the surfaces.
To achieve my goal, I turned to Python (specifically the microgen library of Pr. Yves CHEMISKY). Here is the result (viewed with 3D Viewer) that I would like to obtain using MATLAB:
This result was achieved by performing a Boolean mesh operation between a cylinder and a gyroid. This gave me the idea to check if a similar function existed in MATLAB. I couldn’t find a built-in function, but the mesh_boolean function from the gptoolbox toolbox seems to do the same thing. I tried installing the toolbox, but unfortunately, the function I need does not work, and some libraries are missing. I’ve tried installing them but have not succeeded so far.
Nevertheless, I would like to successfully create this cylindrical gyroid using MATLAB.
Thank you very much for your contribution, and I remain open to other suggestions if anyone has a solution for me.

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!