Draw partial spheroid include a spheroid

I want to draw 1/8 spheroid include a small spheroid and output the geometry for mesh. But my current coding always have discontinue in the cutting plan.

Can anyone help provide a idea of cutting the spheroid in 1/8 not for showing but get the data.

Answer 1

[Bruno Luong]

The code bellow us this FEX to generate mesh points.

% radius of the inner/outer spherical parts
r1 = 1;
r2 = 2;
n = 20; % discretization parameters of spherical parts
W = allVL1(3, n); % FEX file

% Connectivity
XY = W*[0 1 0;
        0 0 1].';
F=delaunay(XY);

% Points in S2
Tri3 = eye(3);
W = W/n;
XYZ = W*Tri3';
XYZ = XYZ ./ sqrt(sum(XYZ.^2,2));

% inner/outer sphericals
XYZ1 = XYZ*r1;
XYZ2 = XYZ*r2;

% TRUE for points on the boundary
ibdr = W==0; 

close all
patcharg = {'FaceColor', 'g', 'FaceAlpha', 0.5};
patch('Faces', F, 'Vertices', XYZ1, patcharg{:});
patch('Faces', F, 'Vertices', XYZ2, patcharg{:});
for k=1:3
    XYZk = [XYZ1(ibdr(:,k),:); flipud(XYZ2(ibdr(:,k),:))];
    % You are free to mesh XYZk, I leave it as polygonal shape (quarter of a rings)
    patch('XData', XYZk(:,1), 'YData', XYZk(:,2), 'ZData', XYZk(:,3), patcharg{:});
end
view(3)
axis equal

[KOU DU]

Thanks Bruno.

In fact, the inner part is a complicated geometry defined by :

abs((x-xc)/a).^(2*p)+abs((y-yc)/b).^(2*p)+abs((z-zc)/c).^(2*p)=1.0

I just show a simple case in my first question. I'll try to change your code to draw this inner part.

One more question. How can I combine all patch results to output?

[Bruno Luong]

So you have some sort of component-scaled L^(2p) ball for the inner face right?

[KOU DU]

yes, the inner face is defined as the equation I show, and the outer face is a sphere.

[Bruno Luong]

Here we go:

% radius of the outer spherical
r1 = 1;
% parameter of inner part
a=0.2;
b=0.2;
c=0.2;
p=0.5;

n = 20; % discretization parameters of spherical parts
W = allVL1(3, n); % FEX file
% Connectivity
XY = W*[0 1 0;
        0 0 1].';
F=delaunay(XY);
% Points in S2
Tri3 = eye(3);
W = W/n;
XYZ = W*Tri3';
XYZ = XYZ ./ sqrt(sum(XYZ.^2,2));
% inner/outer sphericals
XYZ1 = XYZ*r1;

q = 2*p;
qnorm = sum((XYZ ./ [a,b,c]).^q,2).^(1/q); % add abs if components are negative
XYZ2 = XYZ ./ qnorm;
% TRUE for points on the boundary
ibdr = W==0;

% Group vertices
XYZ = [XYZ1; XYZ2];
m = size(XYZ1,1);
Fout = 0 + F;
Fin = m + F ;
Fwall = cell(1,3);
for k=1:3
    bk = find(ibdr(:,k));
    Fwall{k} = [bk(:); m+flip(bk(:))]';
end

close all
patcharg = {'FaceColor', 'g', 'FaceAlpha', 0.5};
patch('Faces', [Fin; Fout], 'Vertices', XYZ, patcharg{:});
patch('Faces', cat(1,Fwall{:}), 'Vertices', XYZ, patcharg{:});
view(3)
axis equal

[KOU DU]

Thanks, Bruno. But there something may be wrong.

The inner surface should not be bigger or smaller when I change p, just the geometry should be changed. That means if we don't change a,b,c, the intersection points with axis should not be changed.

[Bruno Luong]

Ops this is correct norm calculation

pnorm = sum((XYZ ./ [a,b,c]).^(2*p),2).^(1/(2*p));

[KOU DU]

Thank you very much, Bruno.

[KOU DU]

One extra question, could we transfer the patch plan to triangulation form?

[Bruno Luong]

Well I wrote:

"% You are free to mesh XYZk, I leave it as polygonal shape (quarter of a rings)"

You can use meshing tool such as MESH2D on FEX if you want mesh with not too elongated triangles.

Otherwise each three rings can be decomposed in quadrilateral by connecting respective points of the inner/outer boundaries, then each can be split in 2 triangles, but they will be elongated.

I don't know what you want to do with the mesh to decide which one is the best at your place.

[KOU DU]

thank you very much bruno! I want to output all vertices & faces information to stl format. I try use https://www.mathworks.com/matlabcentral/fileexchange/20922-stlwrite-write-ascii-or-binary-stl-files to write stl file. And then use another software to do the mesh. But when I directly output the patch.Faces & patch.Vertices, and import to another software, there always duplicated points erros.

[Bruno Luong]

Well, I'll be kind for once and ended up doing the whole thing for you

 % radius of the outer spherical
r1 = 1;
% parameter of inner part
a=0.2; b=0.2; c=0.2; p=0.7;

n = 20; % discretization parameters of spherical parts
W = allVL1(3, n); % FEX file
% Connectivity
XY = W*[0 1 0;
        0 0 1].';
F = delaunay(XY);
% Points in S2
XYZ = W/n;
% inner/outer sphericals
XYZ1 = XYZ .* (r1./ sqrt(sum(XYZ.^2,2)));
q = 2*p;
qnorm = sum((XYZ ./ [a,b,c]).^q,2).^(1/q); % add abs if components are negative
XYZ2 = XYZ ./ qnorm;
% TRUE for points on the boundary
ibdr = W==0;
% Group vertices
XYZ = [XYZ1; XYZ2];
m = size(XYZ1,1);
Fout = 0 + F;
Fin = m + fliplr(F);
Fwall = cell(1,3);
for k=1:3
    bk = find(ibdr(:,k));
    Fk = [[bk(1:end-1) bk(2:end) bk(1:end-1)]+[0 0 m];
          [bk(2:end) bk(1:end-1) bk(2:end)]+[m m 0]];
    % Re-orienting triangular faces so that they are consistently oriented
    T = reshape(XYZ(Fk,:),[],3,3);        
    N = cross(T(:,2,:)-T(:,1,:),T(:,3,:)-T(:,1,:),3);
    reverse = N(:,:,k) > 0;
    Fk(reverse,:) = fliplr(Fk(reverse,:));
    Fwall{k} = Fk;
end

Faces = cat(1,Fin,Fout,Fwall{:});

close all
patcharg = {'FaceColor', 'g', 'FaceAlpha', 0.7};
patch('Faces', Faces, 'Vertices', XYZ, patcharg{:});
view(3)
axis equal

stlwrite('watermeloon.stl', Faces, XYZ)

[KOU DU]

Thank you!

[Kim]

Since I have a similar problem, I tried to compile this programm but as in my case, I always get the error "Input argument must be a triangulation object."

Could anybody point out, what mistake I've been making?

Answer 2

[darova]

Use patch() to generate planes

clc,clear

R = 10;
r = 3;
t = linspace(0,pi/2,20);
x = [r*cos(t) fliplr(R*cos(t))];
y = [r*sin(t) fliplr(R*sin(t))];
patch(x,y,x*0,'b')
hold on
patch(x,x*0,y,'b')
patch(x*0,x,y,'b')
hold off

alpha(0.5)
view(3)

[KOU DU]

Thanks, darova. But I saw 'patch' is just for showing but can't creat the data to output.

[darova]

What format of data output should be?

[KOU DU]

I just see the patch properties and it can output the coordinates. I'll try to use this function to see whether it works. Thanks.

[KOU DU]

Hello, darova. How can I close the outer sphere surface and inter surface. I try lot of ways using patch but not success.

[KOU DU]

I mean how to draw whole geometry as I show?

[darova]

Use surf() to draw 1/8 of sphere and patch() to create a planes

[KOU DU]

Hello, darova.

in fact, my inner geometry is a complicated geometry and I draw by the following coding.

clc,clear
Radius=1; %Radius of the Outer SPHERE
xc=0; %centre of geometry
yc=0;
zc=0;
a=0.2;
b=0.2;
c=0.2;
p=0.5;
% AxisGrid=(0:1/20:1);
AxisGrid=CSGBOXGRID(1,a,50,50);

f = {@(x,y,z)(abs((x-xc)/Radius).^(2)+abs((y-yc)/Radius).^(2)+abs((z-zc)/Radius).^(2)-1.0^(2))
    @(x,y,z)-(abs((x-xc)/a).^(2*p)+abs((y-yc)/b).^(2*p)+abs((z-zc)/c).^(2*p)-1.0^(2*p))
   };

[x,y,z]=meshgrid(AxisGrid,AxisGrid,AxisGrid);

v1 = f{1}(x,y,z);
v2 = max(v1,f{numel(f)}(x,y,z));
surface=isosurface(x,y,z,v2,0);
surface.facecolor='blue';
count=renderpatch(surface);
grid on
view(7,20)
rotate3d on

I try to add patch to fill the cutting plan but not success. Maybe you can give me some suggestion?

[darova]

Are your shapes (sphere and ellispoid?) always centered? Are their centers in origin?

[KOU DU]

yes

[darova]

Here we go

[darova]

Don't know how introduce 'p' parameter in that form

[KOU DU]

Thanks darova. Sorry I don't see the question of shapes. The inner shape is not only a ellispoid.

[KOU DU]

And as you said. I try to introduce p but not success. Whatever, thank you very much.

Answer 3

[Abhisek Pradhan]

Following code may be used as an alternative to draw a sphere. Theta and Phi can be varied to get the desired result.

R=10;
Phi=linspace(-pi,pi);
Theta=linspace(0,2*pi);
[Phi,Theta]=meshgrid(Phi,Theta);
Z=R*sin(Phi);
X=R*cos(Phi).*cos(Theta);
Y=R*cos(Phi).*sin(Theta);
hSurface  = surf(X,Y,Z);
set(hSurface,'FaceColor',[0 0 1], 'FaceAlpha',0.5,'FaceLighting','gouraud','EdgeColor','none');

Refer meshgrid and surf for more information.

[KOU DU]

Thanks, Pradhan. But I know how to draw a whole sphere or other geometry. The problem I meet now is the discontinue in the cutting plan.