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

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!!!!

Answer 1

[TED MOSBY]

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!

Answer 2

[CHEICK IBRAHIM]

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.