Draw partial spheroid include a spheroid

2019 7 29
3 个回答
回答已采纳
更新时间：2019 10 17
6 次查看（30 天）

0 个投票

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.

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

 采纳的回答

Bruno Luong
Bruno Luong 2019-8-21
编辑：Bruno Luong 2019-8-21

1 个投票

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

13 个评论
显示 11更早的评论 隐藏 11更早的评论

KOU DU
KOU DU 2019-8-22
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
Bruno Luong 2019-8-22
So you have some sort of component-scaled L^(2p) ball for the inner face right?
KOU DU
KOU DU 2019-8-22
yes, the inner face is defined as the equation I show, and the outer face is a sphere.
Bruno Luong
Bruno Luong 2019-8-22
编辑：Bruno Luong 2019-8-22
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
KOU DU 2019-8-22
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
Bruno Luong 2019-8-22
编辑：Bruno Luong 2019-8-22
Ops this is correct norm calculation
pnorm = sum((XYZ ./ [a,b,c]).^(2*p),2).^(1/(2*p));
KOU DU
KOU DU 2019-8-22
Thank you very much, Bruno.
KOU DU
KOU DU 2019-8-22
One extra question, could we transfer the patch plan to triangulation form?
Bruno Luong
Bruno Luong 2019-8-22
编辑：Bruno Luong 2019-8-22
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
KOU DU 2019-8-22
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
Bruno Luong 2019-8-22
编辑：Bruno Luong 2019-8-22
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)
Kim
Kim 2019-10-17
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?

请先登录，再进行评论。

更多回答（2 个）

darova
darova 2019-8-19

1 个投票

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)

13 个评论
显示 11更早的评论 隐藏 11更早的评论

KOU DU
KOU DU 2019-8-19
Thanks, darova. But I saw 'patch' is just for showing but can't creat the data to output.
darova
darova 2019-8-19
What format of data output should be?
KOU DU
KOU DU 2019-8-19
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
KOU DU 2019-8-21
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
KOU DU 2019-8-21
I mean how to draw whole geometry as I show?
darova
darova 2019-8-21
Use surf() to draw 1/8 of sphere and patch() to create a planes
KOU DU
KOU DU 2019-8-21
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
darova 2019-8-21
Are your shapes (sphere and ellispoid?) always centered? Are their centers in origin?
darova
darova 2019-8-22
Don't know how introduce 'p' parameter in that form
KOU DU
KOU DU 2019-8-22
Thanks darova. Sorry I don't see the question of shapes. The inner shape is not only a ellispoid.
KOU DU
KOU DU 2019-8-22
And as you said. I try to introduce p but not success. Whatever, thank you very much.

请先登录，再进行评论。

Abhisek Pradhan
Abhisek Pradhan 2019-8-7

0 个投票

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.

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

KOU DU
KOU DU 2019-8-19
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Surface and Mesh Plots 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by