DEMO_MixedTetHexMeshing
Below is a demonstration of how to create a mixed mesh consisting of (linear) hexahedral and tetrahedral elements. The hexahedral mesh is regular while the tetrahedral mesh is derived using TetGen.
Contents
clear; close all; clc;
Plot settings for the examples below
fontSize=20; faceAlpha1=1; faceAlpha2=0.3; plotColors=gjet(4);
searchRadius=6;
CONVERTING A TRIANGULATED SURFACE TO AN IMAGE WITH DESIRED SIZE, VOXEL SIZE AND ORIGIN
Defining an example triangulated surface model
[Fs,Vs]=stanford_bunny;
Setting control parameters
% Defining the full set of possible control parameters voxelSize=6; % The output image voxel size. imOrigin=min(Vs,[],1)-voxelSize; imMax=max(Vs,[],1)+voxelSize; imSiz=round((imMax-imOrigin)/voxelSize); imSiz=imSiz([2 1 3]); %Image size (x, y corresponds to j,i in image coordinates, hence the permutation) % Using |triSurf2Im| function to convert patch data to image data [M,~]=triSurf2Im(Fs,Vs,voxelSize,imOrigin,imSiz);
Plotting the results
hf1=cFigure; subplot(1,2,1); title('Closed triangulated surface','FontSize',fontSize); xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on; gpatch(Fs,Vs,'g','none'); axis equal; view(3); axis tight; grid on; set(gca,'FontSize',fontSize); camlight('headlight'); lighting phong; set(gca,'fontSize',fontSize); subplot(1,2,2); title('Boundary, intertior and exterior image','FontSize',fontSize); xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on; gpatch(Fs,Vs,'g','none',faceAlpha2); L_plot=false(size(M)); L_plot(:,:,round(size(M,3)/2))=1; [Fm,Vm,Cm]=ind2patch(L_plot,double(M),'sk'); [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3)); Vm=Vm+imOrigin(ones(size(Vm,1),1),:); gpatch(Fm,Vm,Cm,'k'); L_plot=false(size(M));L_plot(round(size(M,1)/2),:,:)=1; [Fm,Vm,Cm]=ind2patch(L_plot,M,'si'); [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3)); Vm=Vm+imOrigin(ones(size(Vm,1),1),:); gpatch(Fm,Vm,Cm,'k'); L_plot=false(size(M));L_plot(:,round(size(M,2)/2),:)=1; [Fm,Vm,Cm]=ind2patch(L_plot,M,'sj'); [Vm(:,1),Vm(:,2),Vm(:,3)]=im2cart(Vm(:,2),Vm(:,1),Vm(:,3),voxelSize*ones(1,3)); Vm=Vm+imOrigin(ones(size(Vm,1),1),:); gpatch(Fm,Vm,Cm,'k'); colormap(gray(3)); caxis([0 2]); hc=colorbar; set(hc,'YTick',[1/3 1 5/3]); set(hc,'YTickLabel',{'Exterior','Boundary','Intertior'}); set(hc,'fontSize',fontSize); axis equal; view(3); axis tight; grid on; set(gca,'FontSize',fontSize); set(gca,'fontSize',fontSize); drawnow;

GET HEXAHEDRAL ELEMENT SET
L_model=(M==2); %Interior&Boundary choosen here %Defining erosion/dilation kernel k=3; p=k-round(k./2); hb=zeros(3,3); hb(2,2,2)=1; hb(2,2,1)=1; hb(2,2,3)=1; hb(1,2,2)=1; hb(3,2,2)=1; hb(2,3,2)=1; hb(2,1,2)=1; L_model_rep=zeros(size(L_model)+(2.*p)); L_model_rep(p+(1:size(L_model,1)),p+(1:size(L_model,2)),p+(1:size(L_model,3)))=L_model; L_model_blur = convn(double(L_model_rep),hb,'valid'); L_model=L_model_blur>=(sum(hb(:))); [E_hex,V_hex,C_hex]=ind2patch(L_model,M,'hu'); % Convert Coordinates [V_hex(:,1),V_hex(:,2),V_hex(:,3)]=im2cart(V_hex(:,2),V_hex(:,1),V_hex(:,3),voxelSize*ones(1,3)); V_hex=V_hex+imOrigin(ones(size(V_hex,1),1),:); % Use element2patch to get patch data to plot the model [F_hex_cut1,C_hex_F]=element2patch(E_hex,C_hex); %Pass through unique_patch to reduce "weight" of plot [Fp,Vp,~,~,~,F_count]=unique_patch(F_hex_cut1,V_hex,[],5); logicUni=F_count==1; %Logic for boundary faces Fq=Fp(logicUni,:); Vq=Vp; [Fq,Vq,~]=patchCleanUnused(Fq,Vq); [Ft,Vt]=quad2tri(Fq,Vq,'b');
Plotting the results
hf1=cFigure; title('Visualizing internal voxels=hexahedral elements','FontSize',fontSize); xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on; gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2); gpatch(Fq,Vq,plotColors(2,:),'k'); camlight('headlight'); lighting flat; axis equal; view(3); axis tight; grid on; set(gca,'FontSize',fontSize); drawnow;

%Joining surface sets F=[Fs;Ft+size(Vs,1)]; V=[Vs;Vt]; C_tet=[ones(size(Fs,1),1);2*ones(size(Ft,1),1)]; %Surface marker colors
Get hole point
[V_hole]=getInnerPoint(Ft,Vt,searchRadius,voxelSize/2,0); plotV(V_hole,'r.','MarkerSize',25);

Get region point
L_in=(M==1); [indInternal]=getInnerVoxel(double(L_in),searchRadius,0); [I_in,J_in,K_in]=ind2sub(size(L_in),indInternal); %Convert to subscript coordinates [X_in,Y_in,Z_in]=im2cart(I_in,J_in,K_in,voxelSize*ones(1,3)); V_in=[X_in Y_in Z_in]; V_in=V_in+imOrigin(ones(size(V_in,1),1),:); plotV(V_in,'k.','MarkerSize',25);

DEFINE FACE BOUNDARY MARKERS
faceBoundaryMarker=C_tet;
Define region points
V_regions=[V_in];
Define hole points
V_holes=[V_hole];
Regional mesh parameters
[edgeLengths]=patchEdgeLengths(F,V);
edgeLengthsMean=mean(edgeLengths);
meanProposedVolume=edgeLengthsMean^3./(6*sqrt(2)); %For regular tetrahedron
regionA=meanProposedVolume;
CREATING THE SMESH STRUCTURE, meshing without the surface constraints imposed by the -Y this time.
stringOpt='-pq1.2AaY'; modelName='tetGenModel'; meshStruct.stringOpt=stringOpt; meshStruct.Faces=F; meshStruct.Nodes=V; meshStruct.holePoints=V_holes; meshStruct.faceBoundaryMarker=faceBoundaryMarker; %Face boundary markers meshStruct.regionPoints=V_regions; %region points meshStruct.regionA=regionA; meshStruct.minRegionMarker=2; %Minimum region marker meshStruct.modelName=modelName;
Mesh model using tetrahedral elements using tetGen (see: http://wias-berlin.de/software/tetgen/)
[meshOutput]=runTetGen(meshStruct); %Run tetGen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- TETGEN Tetrahedral meshing --- 11-Dec-2020 12:15:57 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- Writing SMESH file --- 11-Dec-2020 12:15:57 ----> Adding node field ----> Adding facet field ----> Adding holes specification ----> Adding region specification --- Done --- 11-Dec-2020 12:15:57 --- Running TetGen to mesh input boundary--- 11-Dec-2020 12:15:57 Opening /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.smesh. Delaunizing vertices... Delaunay seconds: 0.041198 Creating surface mesh ... Surface mesh seconds: 0.005787 Recovering boundaries... Boundary recovery seconds: 0.027355 Removing exterior tetrahedra ... Spreading region attributes. Exterior tets removal seconds: 0.006977 Suppressing Steiner points ... Steiner suppression seconds: 1.5e-05 Recovering Delaunayness... Delaunay recovery seconds: 0.011517 Refining mesh... Refinement seconds: 0.039581 Optimizing mesh... Optimization seconds: 0.003861 Jettisoning redundant points. Writing /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.1.node. Writing /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.1.ele. Writing /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.1.face. Writing /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.1.edge. Output seconds: 0.049431 Total running seconds: 0.185938 Statistics: Input points: 2729 Input facets: 5460 Input segments: 8185 Input holes: 1 Input regions: 1 Mesh points: 3514 Mesh tetrahedra: 13888 Mesh faces: 30506 Mesh faces on exterior boundary: 5460 Mesh faces on input facets: 5460 Mesh edges on input segments: 8185 Steiner points inside domain: 785 --- Done --- 11-Dec-2020 12:15:58 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- Importing TetGen files --- 11-Dec-2020 12:15:58 --- Done --- 11-Dec-2020 12:15:58
Access model element and patch data
F_tet_cut1=meshOutput.faces; V_tet=meshOutput.nodes; C_tet=meshOutput.faceMaterialID; E_tet=meshOutput.elements; indBoundary=meshOutput.facesBoundary(meshOutput.boundaryMarker==1,:); indBoundary=unique(indBoundary(:));
MERGING NODE SETS
V=[V_tet;V_hex]; E_hex=E_hex+size(V_tet,1); [~,V,ind1,ind2]=mergeVertices(F_tet_cut1,V); E_tet=ind2(E_tet); E_hex=ind2(E_hex); indBoundary=ind2(indBoundary); E={E_tet, E_hex};
Visualizing mesh
cFigure; title('Mixed TET/HEX mesh','FontSize',fontSize); %Selecting half of the model to see interior X=V(:,2); XE=mean(X(E{1}),2); L=XE>mean(X); [F_tet_cut1,~]=element2patch(E{1}(L,:),C_tet(L)); %Selecting half of the model to see interior X=V(:,2); XE=mean(X(E{2}),2); L=XE>mean(X); [F_hex_cut1,~]=element2patch(E{2}(L,:),C_hex(L)); gpatch(F_tet_cut1,V,plotColors(1,:),'k'); gpatch(F_hex_cut1,V,plotColors(2,:),'k'); gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2); axisGeom(gca,fontSize); camlight headlight; set(gca,'FontSize',fontSize); drawnow;

Smoothing meshes
cPar.Method='LAP';
cPar.n=25;
cPar.RigidConstraints=indBoundary;
[F1,~]=element2patch(E{1},[]);
[F1,~,~]=uniqueIntegerRow(F1);
[F2,~]=element2patch(E{2},[]);
[F2,~,~]=uniqueIntegerRow(F2);
[~,IND_V1,~]=tesIND(F1,V,0);
[~,IND_V2,~]=tesIND(F2,V,0);
IND_V=[IND_V1 IND_V2];
[VS]=tesSmooth([],V,IND_V,cPar);
Visualizing mesh
hf=cFigure; title('Mixed TET/HEX mesh','FontSize',fontSize); %Selecting half of the model to see interior X=V(:,2); XE=mean(X(E{1}),2); L=XE>mean(X); [F_tet_cut1,~]=element2patch(E{1}(L,:),C_tet(L)); %Selecting half of the model to see interior X=V(:,2); XE=mean(X(E{2}),2); L=XE>mean(X); [F_hex_cut1,~]=element2patch(E{2}(L,:),C_hex(L)); hp1=gpatch(F_tet_cut1,VS,plotColors(1,:),'k'); hp2=gpatch(F_hex_cut1,VS,plotColors(2,:),'k'); gpatch(Fs,Vs,0.5*ones(1,3),'none',faceAlpha2); axisGeom(gca,fontSize); camlight headlight; set(gca,'FontSize',fontSize); drawnow;

Set up animation
nSteps=25; %Number of animation steps X=V(:,2); XE1=mean(X(E{1}),2); XE2=mean(X(E{2}),2); animStruct.Time=linspace(0,1,nSteps); %Time vector cutLevel=linspace(min(X(:)),max(X(:)),nSteps); %Property to set for q=1:1:nSteps %Step through time cutLevelNow=cutLevel(q); %The current cut level L1=XE1>cutLevelNow; [F_tet_cut1,~]=element2patch(E{1}(L1,:)); L2=XE2>cutLevelNow; [F_hex_cut1,~]=element2patch(E{2}(L2,:)); %Set entries in animation structure animStruct.Handles{q}=[hp1 hp2]; %Handles of objects to animate animStruct.Props{q}={'Faces','Faces'}; %Properties of objects to animate animStruct.Set{q}={F_tet_cut1,F_hex_cut1}; %Property values for to set in order to animate end %Add animation layer anim8(hf,animStruct);

GIBBON
Kevin M. Moerman ([email protected])
GIBBON footer text
License: https://github.com/gibbonCode/GIBBON/blob/master/LICENSE
GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for image segmentation, image-based modeling, meshing, and finite element analysis.
Copyright (C) 2006-2020 Kevin Mattheus Moerman
This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.