Closed AbrahamKalu-Uka closed 1 year ago
The file above is the stl file I want to mesh with tetGen, Gibbon code
The STL had inconsistent face orientations and was not meshed properly (for finite element analysis), i.e. the surface needs remeshing.
These are the input normals, not the bottom was facing the wrong way:
After reorienting with patchNormalFix
we get:
Now though we see that the top and bottom surfaces have extremely sharp triangles, which is likely not suitable for finite element analysis (and note also that TetGen is not good at remeshing surfaces).
You can remesh the entire surface using something like ggremesh
, however that would not respect your "sharp" edges which I judge are probably very important to your application. So instead you could remesh the top and bottom only for instance, and "hold on" to the boundary nodes to force it to respect your input geometry.
This was the input mesh:
After remeshing we can get something like:
Which is better (and will work with TetGen), but potentially still a bad mesh for finite element analysis.
A similar code can be constructed however to code the mesh entirely in GIBBON and to remesh is fully (top and sides as well as the curves on the boundaries). That would produce an overall homogeneous high quality mesh with a desired point spacing. Let me know if that is what you want. If the side surfaces are simply connecting nodes from the top to the bottom, then this remeshing is easy to do.
Here is the code I used, give it a go:
clear; close all; clc;
%% Control parameters
fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl';
featureAngularTolerance=45/180*pi;
normalDirectionGroupingTolerance=1e-3;
pointSpacing=2;
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces
V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure;
gpatch(F_stl,V_stl); %Visualise surface mesh
patchNormPlot(F_stl,V_stl); %Visualise normal directions
axisGeom; camlight headlight;
gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure;
gpatch(F_stl,V_stl); %Visualise surface mesh
patchNormPlot(F_stl,V_stl); %Visualise normal directions
axisGeom; camlight headlight;
gdrawnow;
%% Feature detect surface
% The following can provide a useful surface labelling. It is repeatable
% for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels
% The top and bottom are found based on the fact that the dot product for
% all face normals should be -1 or 1 with respect to the z-direction
% vector.
numFaceSets=max(C_stl(:));
nz=[0 0 1];
N = patchNormal(F_stl,V_stl);
d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets
dNow=d(C_stl==q);
dDown = abs(dNow+1);
dUp = abs(dNow-1);
if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up
C_top=q;
elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down
C_bottom=q;
end
end
%% Visualise surface labelling
cFigure;
gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh
axisGeom; camlight headlight;
colormap spectral; icolorbar;
gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==C_top,:));
indB_top=edgeListToCurve(Eb_top);
indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==C_bottom,:));
indB_bottom=edgeListToCurve(Eb_bottom);
indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
[F_top,V_top]=regionTriMesh3D({V_stl(indB_top,:)},pointSpacing,0,'natural');
F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_stl(indB_bottom,:)},pointSpacing,0,'natural');
F_bottom=fliplr(F_bottom); %Invert face orientation
%% Visualise new top/bottom meshes
cFigure;
gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
gedge(Eb_top,V_stl,'b',4);
gedge(Eb_bottom,V_stl,'r',4);
gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh
patchNormPlot(F_top,V_top);
gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh
patchNormPlot(F_bottom,V_bottom);
axisGeom; camlight headlight;
gdrawnow;
%% Create a new merged surface model
logicOther=~ismember(C_stl,[C_top, C_bottom]);
F_other=F_stl(logicOther,:);
C_other=C_stl(logicOther);
[F,V,C]=joinElementSets({F_top,F_bottom,F_other},{V_top,V_bottom,V_stl},{C_top.*ones(size(F_top,1),1),C_bottom.*ones(size(F_bottom,1),1),C_other});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
%% Visualise improved mesh
cFigure;
gpatch(F,V,C); %Visualise surface mesh
axisGeom; camlight headlight;
colormap spectral; icolorbar;
gdrawnow;
The code below (I did not have time to comment it very well) does that more parametric approach that I mentioned, and produces an even/homogeneous mesh. See here:
This code first harvests the top and bottom curves from the STL and then proceeds to build geometry from that. An even more reliable code would be to build the geometry (and curves) completely within MATLAB.
Is the milling tool the topic of optimisation? If so, what I suggest is a good idea so you can optimise parameters defining the milling tool design.
clear; close all; clc;
%% Control parameters
fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl';
featureAngularTolerance=45/180*pi;
normalDirectionGroupingTolerance=1e-3;
pointSpacing=0.5;
interpMethod='linear';
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces
V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure;
gpatch(F_stl,V_stl); %Visualise surface mesh
patchNormPlot(F_stl,V_stl); %Visualise normal directions
axisGeom; camlight headlight;
gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure;
gpatch(F_stl,V_stl); %Visualise surface mesh
patchNormPlot(F_stl,V_stl); %Visualise normal directions
axisGeom; camlight headlight;
gdrawnow;
%% Feature detect surface
% The following can provide a useful surface labelling. It is repeatable
% for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels
% The top and bottom are found based on the fact that the dot product for
% all face normals should be -1 or 1 with respect to the z-direction
% vector.
numFaceSets=max(C_stl(:));
nz=[0 0 1];
N = patchNormal(F_stl,V_stl);
d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets
dNow=d(C_stl==q);
dDown = abs(dNow+1);
dUp = abs(dNow-1);
if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up
c_top=q;
elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down
c_bottom=q;
end
end
%% Visualise surface labelling
cFigure;
gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh
axisGeom; camlight headlight;
colormap spectral; icolorbar;
gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==c_top,:));
indB_top=edgeListToCurve(Eb_top);
indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==c_bottom,:));
indB_bottom=edgeListToCurve(Eb_bottom);
indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
logicOther=~ismember(C_stl,[c_top, c_bottom]);
F_other=F_stl(logicOther,:);
C_other=C_stl(logicOther);
c_other=unique(C_other);
%%
CC=patchConnectivity(F_stl,V_stl,'vf');
V2F=CC.vertex.face;
logicValid=V2F>0;
CV=zeros(size(logicValid));
CV(logicValid)=C_stl(V2F(logicValid));
indCornerTop=[];
for q=1:1:numel(indB_top)
cNow=CV(indB_top(q),:);
cNow=unique(cNow(cNow>0));
if any(cNow==c_top) && sum(ismember(cNow,C_other))==2
indCornerTop=[indCornerTop; indB_top(q)];
end
end
indCornerBottom=[];
for q=1:1:numel(indB_bottom)
cNow=CV(indB_bottom(q),:);
cNow=unique(cNow(cNow>0));
if any(cNow==c_bottom) && sum(ismember(cNow,C_other))==2
indCornerBottom=[indCornerBottom; indB_bottom(q)];
end
end
indStartTop=find(indB_top==indCornerTop(1));
if indStartTop>1
indB_top=[indB_top(indStartTop:end) indB_top(1:indStartTop-1)];
end
% Match order
[~,indOrder]=minDist(V_stl(indCornerTop,:),V_stl(indCornerBottom,:));
indCornerBottom=indCornerBottom(indOrder);
indStartBottom=find(indB_bottom==indCornerBottom(1));
if indStartBottom>1
indB_bottom=[indB_bottom(indStartBottom:end) indB_bottom(1:indStartBottom-1)];
end
indOrderCheck=indB_bottom(ismember(indB_bottom,indCornerBottom));
if indOrderCheck(2)~=indCornerBottom(2)
indB_bottom=indB_bottom(numel(indB_bottom):-1:1); %Inverse
indB_bottom=[indB_bottom(end) indB_bottom(1:end-1)];
end
indMust=find(ismember(indB_top,indCornerTop));
V_top_curve=evenlySpaceCurve(V_stl(indB_top,:),pointSpacing,interpMethod,1,indMust);
indMust=find(ismember(indB_bottom,indCornerBottom));
V_bottom_curve=evenlySpaceCurve(V_stl(indB_bottom,:),pointSpacing,interpMethod,1,indMust);
[F_top,V_top]=regionTriMesh3D({V_top_curve},pointSpacing,0,'natural');
F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_bottom_curve},pointSpacing,0,'natural');
% Create coordinate matrices
dm=mean(sqrt(sum((V_top_curve-V_bottom_curve).^2,2)));
numStepsHeight=ceil(dm./pointSpacing);
X=linspacen(V_top_curve(:,1),V_bottom_curve(:,1),numStepsHeight)';
Y=linspacen(V_top_curve(:,2),V_bottom_curve(:,2),numStepsHeight)';
Z=linspacen(V_top_curve(:,3),V_bottom_curve(:,3),numStepsHeight)';
perdiocOpt=[0 1];
[Fqs,Vs]=grid2patch(X,Y,Z,[],perdiocOpt);
Fqs=fliplr(Fqs);
[Fs,Vs]=quad2tri(Fqs,Vs,'f');
%% Visualise
cFigure; hold on;
gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
%plotV(V_stl(indCornerTop,:),'r.','MarkerSize',50)
%plotV(V_stl(indCornerBottom,:),'b.','MarkerSize',50)
plotV(V_top_curve,'r.-','MarkerSize',15,'LineWidth',2)
plotV(V_bottom_curve,'r.-','MarkerSize',15,'LineWidth',2)
% plotV(V_stl(indB_top,:),'r-','LineWidth',2)
% plotV(V_stl(indB_bottom,:),'b-','LineWidth',2)
gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh
patchNormPlot(F_top,V_top);
gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh
patchNormPlot(F_bottom,V_bottom);
gpatch(Fs,Vs,'gw','g',1,1); %Visualise surface mesh
patchNormPlot(Fs,Vs);
axisGeom; camlight headlight;
gdrawnow;
%% Create a new merged surface model
[F,V,C]=joinElementSets({F_top,F_bottom,Fs},{V_top,V_bottom,Vs});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
[Eb]=patchBoundaryLabelEdges(F,V,C);
%% Visualise improved mesh
cFigure;
gpatch(F,V,C); %Visualise surface mesh
gedge(Eb,V,'k',3); %Visualise labelled boundary edges
axisGeom; camlight headlight;
colormap gjet; icolorbar;
gdrawnow;
Hi Kevin,
I have gone through the lines of codes, you sent me and have looked at the screenshot. I am just awestrucked and too excited that I want to speak German. Kevin, du bist Wunderbar, Wunderschön und Super; viel viel viel Dank (you are wonderful, beautiful and great; thank you so so so much). I thank God so much for you and Gibbon code. It was exactly what I needed.
What I need now is what particular information in addition to the nodes (faces or elements) that I need to write into the Abaqus file in order to have a proper tet mesh representation.
Then, I want to carry out optimisation studies considering the endmill shape parameters. And I foresee a challenge. The rake face of the current endmill is at 0 degrees and at the moment I can successfully vary it from 0 to higher negative rake angles without having problems when writing the stl files. However, when I have to vary for positive rake angles, the problem of stl faces overlapping over areas that should not have faces arises. This is because the positive rake face curves inside like a concave surface. And the stl face is written from the outer points (including this concave endpoint) to the center. Please I would like to receive suggestions on how to overcome this challenge.
Like I mentioned in my previous email response to your other mail sent to me, I am away from my office. When I am back I could send a picture showing this problem.
Thank you Kevin and to your Gibbon code colleagues. I am very grateful. I would be sure to reference Gibbon while writing my research papers for publication.
Kindest regards, Abraham Kalu-Uka.
-Ing. Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kevin Mattheus Moerman @.***> Sent: Tuesday, August 1, 2023 4:40:24 PM To: gibbonCode/GIBBON Cc: Kalu-Uka, Abraham; Author Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
The code below (I did not have time to comment it very well) does that more parametric approach that I mentioned, and produces an even/homogeneous mesh. See here: [Screenshot from 2023-08-01 15-36-51]https://user-images.githubusercontent.com/8392709/257549066-50cc53e5-f0dc-48b2-a9ed-226354892bc3.png
This code first harvests the top and bottom curves from the STL and then proceeds to build geometry from that. An even more reliable code would be to build the geometry (and curves) completely within MATLAB.
Is the milling tool the topic of optimisation? If so, what I suggest is a good idea so you can optimise parameters defining the milling tool design.
clear; close all; clc;
%% Control parameters fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl'; featureAngularTolerance=45/180*pi; normalDirectionGroupingTolerance=1e-3; pointSpacing=0.5; interpMethod='linear';
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Feature detect surface % The following can provide a useful surface labelling. It is repeatable % for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels % The top and bottom are found based on the fact that the dot product for % all face normals should be -1 or 1 with respect to the z-direction % vector.
numFaceSets=max(C_stl(:)); nz=[0 0 1]; N = patchNormal(F_stl,V_stl); d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets dNow=d(C_stl==q); dDown = abs(dNow+1); dUp = abs(dNow-1); if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up c_top=q; elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down c_bottom=q; end end
%% Visualise surface labelling
cFigure; gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh axisGeom; camlight headlight; colormap spectral; icolorbar; gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==c_top,:)); indB_top=edgeListToCurve(Eb_top); indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==c_bottom,:)); indB_bottom=edgeListToCurve(Eb_bottom); indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
logicOther=~ismember(C_stl,[c_top, c_bottom]); F_other=F_stl(logicOther,:); C_other=C_stl(logicOther); c_other=unique(C_other);
%%
CC=patchConnectivity(F_stl,V_stl,'vf');
V2F=CC.vertex.face; logicValid=V2F>0; CV=zeros(size(logicValid)); CV(logicValid)=C_stl(V2F(logicValid));
indCornerTop=[]; for q=1:1:numel(indB_top) cNow=CV(indB_top(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_top) && sum(ismember(cNow,C_other))==2 indCornerTop=[indCornerTop; indB_top(q)]; end end
indCornerBottom=[]; for q=1:1:numel(indB_bottom) cNow=CV(indB_bottom(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_bottom) && sum(ismember(cNow,C_other))==2 indCornerBottom=[indCornerBottom; indB_bottom(q)]; end end
indStartTop=find(indB_top==indCornerTop(1)); if indStartTop>1 indB_top=[indB_top(indStartTop:end) indB_top(1:indStartTop-1)]; end
% Match order [~,indOrder]=minDist(V_stl(indCornerTop,:),V_stl(indCornerBottom,:)); indCornerBottom=indCornerBottom(indOrder);
indStartBottom=find(indB_bottom==indCornerBottom(1)); if indStartBottom>1 indB_bottom=[indB_bottom(indStartBottom:end) indB_bottom(1:indStartBottom-1)]; end
indOrderCheck=indB_bottom(ismember(indB_bottom,indCornerBottom)); if indOrderCheck(2)~=indCornerBottom(2) indB_bottom=indB_bottom(numel(indB_bottom):-1:1); %Inverse indB_bottom=[indB_bottom(end) indB_bottom(1:end-1)]; end
indMust=find(ismember(indB_top,indCornerTop)); V_top_curve=evenlySpaceCurve(V_stl(indB_top,:),pointSpacing,interpMethod,1,indMust);
indMust=find(ismember(indB_bottom,indCornerBottom)); V_bottom_curve=evenlySpaceCurve(V_stl(indB_bottom,:),pointSpacing,interpMethod,1,indMust);
[F_top,V_top]=regionTriMesh3D({V_top_curve},pointSpacing,0,'natural'); F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_bottom_curve},pointSpacing,0,'natural');
% Create coordinate matrices dm=mean(sqrt(sum((V_top_curve-V_bottom_curve).^2,2))); numStepsHeight=ceil(dm./pointSpacing);
X=linspacen(V_top_curve(:,1),V_bottom_curve(:,1),numStepsHeight)'; Y=linspacen(V_top_curve(:,2),V_bottom_curve(:,2),numStepsHeight)'; Z=linspacen(V_top_curve(:,3),V_bottom_curve(:,3),numStepsHeight)';
perdiocOpt=[0 1]; [Fqs,Vs]=grid2patch(X,Y,Z,[],perdiocOpt); Fqs=fliplr(Fqs); [Fs,Vs]=quad2tri(Fqs,Vs,'f');
%% Visualise
cFigure; hold on; gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
%plotV(V_stl(indCornerTop,:),'r.','MarkerSize',50) %plotV(V_stl(indCornerBottom,:),'b.','MarkerSize',50) plotV(V_top_curve,'r.-','MarkerSize',15,'LineWidth',2) plotV(V_bottom_curve,'r.-','MarkerSize',15,'LineWidth',2)
% plotV(V_stl(indB_top,:),'r-','LineWidth',2) % plotV(V_stl(indB_bottom,:),'b-','LineWidth',2) gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh patchNormPlot(F_top,V_top); gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh patchNormPlot(F_bottom,V_bottom); gpatch(Fs,Vs,'gw','g',1,1); %Visualise surface mesh patchNormPlot(Fs,Vs); axisGeom; camlight headlight; gdrawnow;
%% Create a new merged surface model
[F,V,C]=joinElementSets({F_top,F_bottom,Fs},{V_top,V_bottom,Vs});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
[Eb]=patchBoundaryLabelEdges(F,V,C);
%% Visualise improved mesh
cFigure; gpatch(F,V,C); %Visualise surface mesh gedge(Eb,V,'k',3); %Visualise labelled boundary edges
axisGeom; camlight headlight; colormap gjet; icolorbar; gdrawnow;
— Reply to this email directly, view it on GitHubhttps://github.com/gibbonCode/GIBBON/issues/170#issuecomment-1660465560, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BADVYKDWQ6FUSRHBRJWE26LXTEINRANCNFSM6AAAAAAZEQTQY4. You are receiving this because you authored the thread.Message ID: @.***>
Good day Kevin,
I am now back to the office and in this email, I have attached stl and Abaqus input files showing the different
challenges I have to face while using Gibbon Code.
I look forward to receiving guidance and suggestions on how to solve this problem.
Kind regards, Abraham Kalu-Uka.
-Ing. Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/http://www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kalu-Uka, Abraham Sent: Wednesday, August 2, 2023 10:06:24 PM To: gibbonCode/GIBBON Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
Hi Kevin,
I have gone through the lines of codes, you sent me and have looked at the screenshot. I am just awestrucked and too excited that I want to speak German. Kevin, du bist Wunderbar, Wunderschön und Super; viel viel viel Dank (you are wonderful, beautiful and great; thank you so so so much). I thank God so much for you and Gibbon code. It was exactly what I needed.
What I need now is what particular information in addition to the nodes (faces or elements) that I need to write into the Abaqus file in order to have a proper tet mesh representation.
Then, I want to carry out optimisation studies considering the endmill shape parameters. And I foresee a challenge. The rake face of the current endmill is at 0 degrees and at the moment I can successfully vary it from 0 to higher negative rake angles without having problems when writing the stl files. However, when I have to vary for positive rake angles, the problem of stl faces overlapping over areas that should not have faces arises. This is because the positive rake face curves inside like a concave surface. And the stl face is written from the outer points (including this concave endpoint) to the center. Please I would like to receive suggestions on how to overcome this challenge.
Like I mentioned in my previous email response to your other mail sent to me, I am away from my office. When I am back I could send a picture showing this problem.
Thank you Kevin and to your Gibbon code colleagues. I am very grateful. I would be sure to reference Gibbon while writing my research papers for publication.
Kindest regards, Abraham Kalu-Uka.
-Ing. Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kevin Mattheus Moerman @.***> Sent: Tuesday, August 1, 2023 4:40:24 PM To: gibbonCode/GIBBON Cc: Kalu-Uka, Abraham; Author Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
The code below (I did not have time to comment it very well) does that more parametric approach that I mentioned, and produces an even/homogeneous mesh. See here: [Screenshot from 2023-08-01 15-36-51]https://user-images.githubusercontent.com/8392709/257549066-50cc53e5-f0dc-48b2-a9ed-226354892bc3.png
This code first harvests the top and bottom curves from the STL and then proceeds to build geometry from that. An even more reliable code would be to build the geometry (and curves) completely within MATLAB.
Is the milling tool the topic of optimisation? If so, what I suggest is a good idea so you can optimise parameters defining the milling tool design.
clear; close all; clc;
%% Control parameters fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl'; featureAngularTolerance=45/180*pi; normalDirectionGroupingTolerance=1e-3; pointSpacing=0.5; interpMethod='linear';
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Feature detect surface % The following can provide a useful surface labelling. It is repeatable % for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels % The top and bottom are found based on the fact that the dot product for % all face normals should be -1 or 1 with respect to the z-direction % vector.
numFaceSets=max(C_stl(:)); nz=[0 0 1]; N = patchNormal(F_stl,V_stl); d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets dNow=d(C_stl==q); dDown = abs(dNow+1); dUp = abs(dNow-1); if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up c_top=q; elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down c_bottom=q; end end
%% Visualise surface labelling
cFigure; gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh axisGeom; camlight headlight; colormap spectral; icolorbar; gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==c_top,:)); indB_top=edgeListToCurve(Eb_top); indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==c_bottom,:)); indB_bottom=edgeListToCurve(Eb_bottom); indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
logicOther=~ismember(C_stl,[c_top, c_bottom]); F_other=F_stl(logicOther,:); C_other=C_stl(logicOther); c_other=unique(C_other);
%%
CC=patchConnectivity(F_stl,V_stl,'vf');
V2F=CC.vertex.face; logicValid=V2F>0; CV=zeros(size(logicValid)); CV(logicValid)=C_stl(V2F(logicValid));
indCornerTop=[]; for q=1:1:numel(indB_top) cNow=CV(indB_top(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_top) && sum(ismember(cNow,C_other))==2 indCornerTop=[indCornerTop; indB_top(q)]; end end
indCornerBottom=[]; for q=1:1:numel(indB_bottom) cNow=CV(indB_bottom(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_bottom) && sum(ismember(cNow,C_other))==2 indCornerBottom=[indCornerBottom; indB_bottom(q)]; end end
indStartTop=find(indB_top==indCornerTop(1)); if indStartTop>1 indB_top=[indB_top(indStartTop:end) indB_top(1:indStartTop-1)]; end
% Match order [~,indOrder]=minDist(V_stl(indCornerTop,:),V_stl(indCornerBottom,:)); indCornerBottom=indCornerBottom(indOrder);
indStartBottom=find(indB_bottom==indCornerBottom(1)); if indStartBottom>1 indB_bottom=[indB_bottom(indStartBottom:end) indB_bottom(1:indStartBottom-1)]; end
indOrderCheck=indB_bottom(ismember(indB_bottom,indCornerBottom)); if indOrderCheck(2)~=indCornerBottom(2) indB_bottom=indB_bottom(numel(indB_bottom):-1:1); %Inverse indB_bottom=[indB_bottom(end) indB_bottom(1:end-1)]; end
indMust=find(ismember(indB_top,indCornerTop)); V_top_curve=evenlySpaceCurve(V_stl(indB_top,:),pointSpacing,interpMethod,1,indMust);
indMust=find(ismember(indB_bottom,indCornerBottom)); V_bottom_curve=evenlySpaceCurve(V_stl(indB_bottom,:),pointSpacing,interpMethod,1,indMust);
[F_top,V_top]=regionTriMesh3D({V_top_curve},pointSpacing,0,'natural'); F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_bottom_curve},pointSpacing,0,'natural');
% Create coordinate matrices dm=mean(sqrt(sum((V_top_curve-V_bottom_curve).^2,2))); numStepsHeight=ceil(dm./pointSpacing);
X=linspacen(V_top_curve(:,1),V_bottom_curve(:,1),numStepsHeight)'; Y=linspacen(V_top_curve(:,2),V_bottom_curve(:,2),numStepsHeight)'; Z=linspacen(V_top_curve(:,3),V_bottom_curve(:,3),numStepsHeight)';
perdiocOpt=[0 1]; [Fqs,Vs]=grid2patch(X,Y,Z,[],perdiocOpt); Fqs=fliplr(Fqs); [Fs,Vs]=quad2tri(Fqs,Vs,'f');
%% Visualise
cFigure; hold on; gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
%plotV(V_stl(indCornerTop,:),'r.','MarkerSize',50) %plotV(V_stl(indCornerBottom,:),'b.','MarkerSize',50) plotV(V_top_curve,'r.-','MarkerSize',15,'LineWidth',2) plotV(V_bottom_curve,'r.-','MarkerSize',15,'LineWidth',2)
% plotV(V_stl(indB_top,:),'r-','LineWidth',2) % plotV(V_stl(indB_bottom,:),'b-','LineWidth',2) gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh patchNormPlot(F_top,V_top); gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh patchNormPlot(F_bottom,V_bottom); gpatch(Fs,Vs,'gw','g',1,1); %Visualise surface mesh patchNormPlot(Fs,Vs); axisGeom; camlight headlight; gdrawnow;
%% Create a new merged surface model
[F,V,C]=joinElementSets({F_top,F_bottom,Fs},{V_top,V_bottom,Vs});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
[Eb]=patchBoundaryLabelEdges(F,V,C);
%% Visualise improved mesh
cFigure; gpatch(F,V,C); %Visualise surface mesh gedge(Eb,V,'k',3); %Visualise labelled boundary edges
axisGeom; camlight headlight; colormap gjet; icolorbar; gdrawnow;
— Reply to this email directly, view it on GitHubhttps://github.com/gibbonCode/GIBBON/issues/170#issuecomment-1660465560, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BADVYKDWQ6FUSRHBRJWE26LXTEINRANCNFSM6AAAAAAZEQTQY4. You are receiving this because you authored the thread.Message ID: @.***>
Dear Kevin,
The milling tool is the topic for optimization. My plan is to use this end milling tool to simulate a milling process using Abaqus. I get the
average cutting force and focus on minimizing this cutting force. Originally, I used a Bspline to generate the profile of the endmill tool after
which I wrote an stl file for the faces and vertices. The Bspline used points which were generated from analytical equations that describe the
various angles that characterize an endmill.
My present difficulty is establishing the constraints (e.g. the curve always being a monotonic curve) on the endmill curve in order to have a smooth meaningful curve each time the parametric
variables are changed. I also want to keep the mesh details constant especially the node that is defined as the Reference Point for turning the endmill.
Please I would greatly appreciate whatever guidance you would give me.
Kind regards,
Abraham.
Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/http://www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kevin Mattheus Moerman @.***> Sent: Tuesday, August 1, 2023 4:40:24 PM To: gibbonCode/GIBBON Cc: Kalu-Uka, Abraham; Author Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
The code below (I did not have time to comment it very well) does that more parametric approach that I mentioned, and produces an even/homogeneous mesh. See here: [Screenshot from 2023-08-01 15-36-51]https://user-images.githubusercontent.com/8392709/257549066-50cc53e5-f0dc-48b2-a9ed-226354892bc3.png
This code first harvests the top and bottom curves from the STL and then proceeds to build geometry from that. An even more reliable code would be to build the geometry (and curves) completely within MATLAB.
Is the milling tool the topic of optimisation? If so, what I suggest is a good idea so you can optimise parameters defining the milling tool design.
clear; close all; clc;
%% Control parameters fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl'; featureAngularTolerance=45/180*pi; normalDirectionGroupingTolerance=1e-3; pointSpacing=0.5; interpMethod='linear';
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Feature detect surface % The following can provide a useful surface labelling. It is repeatable % for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels % The top and bottom are found based on the fact that the dot product for % all face normals should be -1 or 1 with respect to the z-direction % vector.
numFaceSets=max(C_stl(:)); nz=[0 0 1]; N = patchNormal(F_stl,V_stl); d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets dNow=d(C_stl==q); dDown = abs(dNow+1); dUp = abs(dNow-1); if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up c_top=q; elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down c_bottom=q; end end
%% Visualise surface labelling
cFigure; gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh axisGeom; camlight headlight; colormap spectral; icolorbar; gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==c_top,:)); indB_top=edgeListToCurve(Eb_top); indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==c_bottom,:)); indB_bottom=edgeListToCurve(Eb_bottom); indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
logicOther=~ismember(C_stl,[c_top, c_bottom]); F_other=F_stl(logicOther,:); C_other=C_stl(logicOther); c_other=unique(C_other);
%%
CC=patchConnectivity(F_stl,V_stl,'vf');
V2F=CC.vertex.face; logicValid=V2F>0; CV=zeros(size(logicValid)); CV(logicValid)=C_stl(V2F(logicValid));
indCornerTop=[]; for q=1:1:numel(indB_top) cNow=CV(indB_top(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_top) && sum(ismember(cNow,C_other))==2 indCornerTop=[indCornerTop; indB_top(q)]; end end
indCornerBottom=[]; for q=1:1:numel(indB_bottom) cNow=CV(indB_bottom(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_bottom) && sum(ismember(cNow,C_other))==2 indCornerBottom=[indCornerBottom; indB_bottom(q)]; end end
indStartTop=find(indB_top==indCornerTop(1)); if indStartTop>1 indB_top=[indB_top(indStartTop:end) indB_top(1:indStartTop-1)]; end
% Match order [~,indOrder]=minDist(V_stl(indCornerTop,:),V_stl(indCornerBottom,:)); indCornerBottom=indCornerBottom(indOrder);
indStartBottom=find(indB_bottom==indCornerBottom(1)); if indStartBottom>1 indB_bottom=[indB_bottom(indStartBottom:end) indB_bottom(1:indStartBottom-1)]; end
indOrderCheck=indB_bottom(ismember(indB_bottom,indCornerBottom)); if indOrderCheck(2)~=indCornerBottom(2) indB_bottom=indB_bottom(numel(indB_bottom):-1:1); %Inverse indB_bottom=[indB_bottom(end) indB_bottom(1:end-1)]; end
indMust=find(ismember(indB_top,indCornerTop)); V_top_curve=evenlySpaceCurve(V_stl(indB_top,:),pointSpacing,interpMethod,1,indMust);
indMust=find(ismember(indB_bottom,indCornerBottom)); V_bottom_curve=evenlySpaceCurve(V_stl(indB_bottom,:),pointSpacing,interpMethod,1,indMust);
[F_top,V_top]=regionTriMesh3D({V_top_curve},pointSpacing,0,'natural'); F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_bottom_curve},pointSpacing,0,'natural');
% Create coordinate matrices dm=mean(sqrt(sum((V_top_curve-V_bottom_curve).^2,2))); numStepsHeight=ceil(dm./pointSpacing);
X=linspacen(V_top_curve(:,1),V_bottom_curve(:,1),numStepsHeight)'; Y=linspacen(V_top_curve(:,2),V_bottom_curve(:,2),numStepsHeight)'; Z=linspacen(V_top_curve(:,3),V_bottom_curve(:,3),numStepsHeight)';
perdiocOpt=[0 1]; [Fqs,Vs]=grid2patch(X,Y,Z,[],perdiocOpt); Fqs=fliplr(Fqs); [Fs,Vs]=quad2tri(Fqs,Vs,'f');
%% Visualise
cFigure; hold on; gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
%plotV(V_stl(indCornerTop,:),'r.','MarkerSize',50) %plotV(V_stl(indCornerBottom,:),'b.','MarkerSize',50) plotV(V_top_curve,'r.-','MarkerSize',15,'LineWidth',2) plotV(V_bottom_curve,'r.-','MarkerSize',15,'LineWidth',2)
% plotV(V_stl(indB_top,:),'r-','LineWidth',2) % plotV(V_stl(indB_bottom,:),'b-','LineWidth',2) gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh patchNormPlot(F_top,V_top); gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh patchNormPlot(F_bottom,V_bottom); gpatch(Fs,Vs,'gw','g',1,1); %Visualise surface mesh patchNormPlot(Fs,Vs); axisGeom; camlight headlight; gdrawnow;
%% Create a new merged surface model
[F,V,C]=joinElementSets({F_top,F_bottom,Fs},{V_top,V_bottom,Vs});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
[Eb]=patchBoundaryLabelEdges(F,V,C);
%% Visualise improved mesh
cFigure; gpatch(F,V,C); %Visualise surface mesh gedge(Eb,V,'k',3); %Visualise labelled boundary edges
axisGeom; camlight headlight; colormap gjet; icolorbar; gdrawnow;
— Reply to this email directly, view it on GitHubhttps://github.com/gibbonCode/GIBBON/issues/170#issuecomment-1660465560, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BADVYKDWQ6FUSRHBRJWE26LXTEINRANCNFSM6AAAAAAZEQTQY4. You are receiving this because you authored the thread.Message ID: @.***>
Dear Kevin,
Please I have a question in addition to the previous inquiries I recently made (even though I have not received
any response yet on them from you).
In the creation of the mesh of the end mill, i want to vary the spline curve that bring about the shape. However, I am
afraid the mesh information would change with each change in the spline curve during an optimization process.
My questions are: 1. Is there any way to control the mesh information to be constant regardless of the change of shape. Please
it is important to point out that the inner and outer radius as well as the height of the tool stays thesame.
flat ends would have thesame node number and position regardless of the change in the mesh during optimization. This is because
I use only one of the 2 center points (node) of a flat endmill as my Reference point during Abaqus simulation. And this node which the
reference point is used to fix the rotation boundary conditions for the endmill. Thus, it has to be constant (in node number and position) so that during simulation,
Abaqus refers to thesame node as the reference point.
Thank you as I look forward to receiving your valuable advice and solution.
Kind regards,
Abraham Kalu-Uka.
Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/http://www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kalu-Uka, Abraham Sent: Friday, October 6, 2023 3:21:45 PM To: gibbonCode/GIBBON Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
Dear Kevin,
The milling tool is the topic for optimization. My plan is to use this end milling tool to simulate a milling process using Abaqus. I get the
average cutting force and focus on minimizing this cutting force. Originally, I used a Bspline to generate the profile of the endmill tool after
which I wrote an stl file for the faces and vertices. The Bspline used points which were generated from analytical equations that describe the
various angles that characterize an endmill.
My present difficulty is establishing the constraints (e.g. the curve always being a monotonic curve) on the endmill curve in order to have a smooth meaningful curve each time the parametric
variables are changed. I also want to keep the mesh details constant especially the node that is defined as the Reference Point for turning the endmill.
Please I would greatly appreciate whatever guidance you would give me.
Kind regards,
Abraham.
Abraham Kalu-Uka, M.Sc. Institute of Engineering and Computational Mechanics (name in German: Institut für Technische und Numerische Mechanik) University of Stuttgart, Pfaffenwaldring 9, 70569 Stuttgart, Germany @.***, tel +49 711-685-64724, fax -66400 www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/http://www.itm.uni-stuttgart.de/institut/team/Kalu-Uka/
From: Kevin Mattheus Moerman @.***> Sent: Tuesday, August 1, 2023 4:40:24 PM To: gibbonCode/GIBBON Cc: Kalu-Uka, Abraham; Author Subject: Re: [gibbonCode/GIBBON] Errors from using the Gibbon Code (Issue #170)
The code below (I did not have time to comment it very well) does that more parametric approach that I mentioned, and produces an even/homogeneous mesh. See here: [Screenshot from 2023-08-01 15-36-51]https://user-images.githubusercontent.com/8392709/257549066-50cc53e5-f0dc-48b2-a9ed-226354892bc3.png
This code first harvests the top and bottom curves from the STL and then proceeds to build geometry from that. An even more reliable code would be to build the geometry (and curves) completely within MATLAB.
Is the milling tool the topic of optimisation? If so, what I suggest is a good idea so you can optimise parameters defining the milling tool design.
clear; close all; clc;
%% Control parameters fileNameSTL='/home/kevin/Downloads/StlMeshProblem(1)/Endmill.stl'; featureAngularTolerance=45/180*pi; normalDirectionGroupingTolerance=1e-3; pointSpacing=0.5; interpMethod='linear';
%% Import surface model geometry
TR = stlread(fileNameSTL); % Read to MATLAB triangulation
F_stl=TR.ConnectivityList; % Get faces V_stl=TR.Points; % Get vertices
[F_stl,V_stl]=mergeVertices(F_stl,V_stl); % Merge shared nodes (required for STL)
% pointSpacing=mean(patchEdgeLengths(F_stl,V_stl)); %Derived control parameter
%% Figure to show that the model has inconsistent normal directions (see bottom face)
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Fix normal directions
[F_stl]=patchNormalFix(F_stl);
%% Figure to show that the model now has consistent normal directions
cFigure; gpatch(F_stl,V_stl); %Visualise surface mesh patchNormPlot(F_stl,V_stl); %Visualise normal directions axisGeom; camlight headlight; gdrawnow;
%% Feature detect surface % The following can provide a useful surface labelling. It is repeatable % for a single model, but labelling may depend on the model used.
[C_stl]=patchFeatureDetect(F_stl,V_stl,featureAngularTolerance);
%% Discover automatically, rather than hardcoded/manually, the top and bottom labels % The top and bottom are found based on the fact that the dot product for % all face normals should be -1 or 1 with respect to the z-direction % vector.
numFaceSets=max(C_stl(:)); nz=[0 0 1]; N = patchNormal(F_stl,V_stl); d=dot(N,nz(ones(size(N,1),1),:),2); %Dot products of normals with Z-direction
for q=1:1:numFaceSets dNow=d(C_stl==q); dDown = abs(dNow+1); dUp = abs(dNow-1); if all(dUp<normalDirectionGroupingTolerance) % All were close to 1, hence up c_top=q; elseif all(dDown<normalDirectionGroupingTolerance) % All were close to -1, hence down c_bottom=q; end end
%% Visualise surface labelling
cFigure; gpatch(F_stl,V_stl,C_stl); %Visualise surface mesh axisGeom; camlight headlight; colormap spectral; icolorbar; gdrawnow;
%% Replace top and bottom surfaces with improved mesh
Eb_top=patchBoundary(F_stl(C_stl==c_top,:)); indB_top=edgeListToCurve(Eb_top); indB_top=indB_top(1:end-1); %end=start hence remove double from list
Eb_bottom=patchBoundary(F_stl(C_stl==c_bottom,:)); indB_bottom=edgeListToCurve(Eb_bottom); indB_bottom=indB_bottom(1:end-1); %end=start hence remove double from list
logicOther=~ismember(C_stl,[c_top, c_bottom]); F_other=F_stl(logicOther,:); C_other=C_stl(logicOther); c_other=unique(C_other);
%%
CC=patchConnectivity(F_stl,V_stl,'vf');
V2F=CC.vertex.face; logicValid=V2F>0; CV=zeros(size(logicValid)); CV(logicValid)=C_stl(V2F(logicValid));
indCornerTop=[]; for q=1:1:numel(indB_top) cNow=CV(indB_top(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_top) && sum(ismember(cNow,C_other))==2 indCornerTop=[indCornerTop; indB_top(q)]; end end
indCornerBottom=[]; for q=1:1:numel(indB_bottom) cNow=CV(indB_bottom(q),:); cNow=unique(cNow(cNow>0)); if any(cNow==c_bottom) && sum(ismember(cNow,C_other))==2 indCornerBottom=[indCornerBottom; indB_bottom(q)]; end end
indStartTop=find(indB_top==indCornerTop(1)); if indStartTop>1 indB_top=[indB_top(indStartTop:end) indB_top(1:indStartTop-1)]; end
% Match order [~,indOrder]=minDist(V_stl(indCornerTop,:),V_stl(indCornerBottom,:)); indCornerBottom=indCornerBottom(indOrder);
indStartBottom=find(indB_bottom==indCornerBottom(1)); if indStartBottom>1 indB_bottom=[indB_bottom(indStartBottom:end) indB_bottom(1:indStartBottom-1)]; end
indOrderCheck=indB_bottom(ismember(indB_bottom,indCornerBottom)); if indOrderCheck(2)~=indCornerBottom(2) indB_bottom=indB_bottom(numel(indB_bottom):-1:1); %Inverse indB_bottom=[indB_bottom(end) indB_bottom(1:end-1)]; end
indMust=find(ismember(indB_top,indCornerTop)); V_top_curve=evenlySpaceCurve(V_stl(indB_top,:),pointSpacing,interpMethod,1,indMust);
indMust=find(ismember(indB_bottom,indCornerBottom)); V_bottom_curve=evenlySpaceCurve(V_stl(indB_bottom,:),pointSpacing,interpMethod,1,indMust);
[F_top,V_top]=regionTriMesh3D({V_top_curve},pointSpacing,0,'natural'); F_top=fliplr(F_top); %Invert face orientation
[F_bottom,V_bottom]=regionTriMesh3D({V_bottom_curve},pointSpacing,0,'natural');
% Create coordinate matrices dm=mean(sqrt(sum((V_top_curve-V_bottom_curve).^2,2))); numStepsHeight=ceil(dm./pointSpacing);
X=linspacen(V_top_curve(:,1),V_bottom_curve(:,1),numStepsHeight)'; Y=linspacen(V_top_curve(:,2),V_bottom_curve(:,2),numStepsHeight)'; Z=linspacen(V_top_curve(:,3),V_bottom_curve(:,3),numStepsHeight)';
perdiocOpt=[0 1]; [Fqs,Vs]=grid2patch(X,Y,Z,[],perdiocOpt); Fqs=fliplr(Fqs); [Fs,Vs]=quad2tri(Fqs,Vs,'f');
%% Visualise
cFigure; hold on; gpatch(F_stl,V_stl,'w','none',0.5); %Visualise surface mesh
%plotV(V_stl(indCornerTop,:),'r.','MarkerSize',50) %plotV(V_stl(indCornerBottom,:),'b.','MarkerSize',50) plotV(V_top_curve,'r.-','MarkerSize',15,'LineWidth',2) plotV(V_bottom_curve,'r.-','MarkerSize',15,'LineWidth',2)
% plotV(V_stl(indB_top,:),'r-','LineWidth',2) % plotV(V_stl(indB_bottom,:),'b-','LineWidth',2) gpatch(F_top,V_top,'b','b',1,2); %Visualise surface mesh patchNormPlot(F_top,V_top); gpatch(F_bottom,V_bottom,'r','r',1,2); %Visualise surface mesh patchNormPlot(F_bottom,V_bottom); gpatch(Fs,Vs,'gw','g',1,1); %Visualise surface mesh patchNormPlot(Fs,Vs); axisGeom; camlight headlight; gdrawnow;
%% Create a new merged surface model
[F,V,C]=joinElementSets({F_top,F_bottom,Fs},{V_top,V_bottom,Vs});
[F,V]=mergeVertices(F,V);
[F,V]=patchCleanUnused(F,V);
[Eb]=patchBoundaryLabelEdges(F,V,C);
%% Visualise improved mesh
cFigure; gpatch(F,V,C); %Visualise surface mesh gedge(Eb,V,'k',3); %Visualise labelled boundary edges
axisGeom; camlight headlight; colormap gjet; icolorbar; gdrawnow;
— Reply to this email directly, view it on GitHubhttps://github.com/gibbonCode/GIBBON/issues/170#issuecomment-1660465560, or unsubscribehttps://github.com/notifications/unsubscribe-auth/BADVYKDWQ6FUSRHBRJWE26LXTEINRANCNFSM6AAAAAAZEQTQY4. You are receiving this because you authored the thread.Message ID: @.***>
We dealt with this in a video call, so I will close this for now.
Hello Kevin,
Thanks for your contributions towards creating the Gibbon code. I have tried using it but I am getting some errors. I did put in alot of effort to find where the problem is, without success. The errors are Index exceeds the number of array elements. Index must not exceed 0.
Error in interp1 (line 153) extrapMask = Xq < X(1) | Xq > X(end);
Error in resampleColormap (line 43) cMap_i(:,q)=interp1(ind,cMap(:,q),ind_i,'linear');
Error in icolorbar (line 74) cn=resampleColormap(c,numel(hc.Ticks));
Error in meshView (line 129) icolorbar;
Error in close_curve (line 468) meshView(meshOutput,[]);
Actually I used the following lines of code (below) to generate a tetmesh from an stl file, and there is a valid inputstruct, but the meshoutput is empty. inputStruct.stringOpt='-pq1.2AaYc'; inputStruct.Faces=fliplr(F); inputStruct.Nodes=V; inputStruct.holePoints=[]; %holes inputStruct.faceBoundaryMarker=ones(size(F,1),1); %Face boundary markers inputStruct.regionPoints=getInnerPoint(F,V); %region interior point inputStruct.regionA=tetVolMeanEst(F,V); %Volume attribute inputStruct.minRegionMarker=2; %Minimum region marker
% Mesh model using tetrahedral elements using tetGen [meshOutput]=runTetGen(inputStruct); %Run tetGen %% meshView(meshOutput,[]);
I would be grateful for your advice and suggestions on how to solve this problem.
Best regards, Abraham.