# DEMO_MixedTetHexMeshing_torus

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;
edgeColor=0.3*ones(1,3);
edgeWidth=1;
```
```searchRadius=6;
```

## CONVERTING A TRIANGULATED SURFACE TO AN IMAGE WITH DESIRED SIZE, VOXEL SIZE AND ORIGIN

Defining an example triangulated surface model

```% Defining a deformed and rotated torus shape
nr=25;
nc=50;
ptype='tri';
[Fs,Vs]=patchTorus(r,nr,rc,nc,ptype);
[THETA,RHO] = cart2pol(Vs(:,1),Vs(:,2));
Vs(:,3)=Vs(:,3)+sin(3*THETA);
[R,~]=euler2DCM([0.5*pi 0.5*pi 0.*pi]);
Vs=Vs*R;
```

Setting control parameters

```% Defining the full set of possible control parameters
voxelSize=0.2; % 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;
patch('Faces',Fs,'Vertices',Vs,'FaceColor','g','EdgeColor','k','FaceAlpha',faceAlpha1);
axis equal; view(3); axis tight;  grid on;  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;

patch('Faces',Fs,'Vertices',Vs,'FaceColor','g','EdgeColor','none','FaceAlpha',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),:);
patch('Faces',Fm,'Vertices',Vm,'FaceColor','flat','CData',Cm,'EdgeColor','k','FaceAlpha',faceAlpha1);

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),:);
patch('Faces',Fm,'Vertices',Vm,'FaceColor','flat','CData',Cm,'EdgeColor','k','FaceAlpha',faceAlpha1);

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),:);
patch('Faces',Fm,'Vertices',Vm,'FaceColor','flat','CData',Cm,'EdgeColor','k','FaceAlpha',faceAlpha1);

colormap(gray(3)); caxis([0 2]);
hc=colorbar;
set(hc,'YTick',[1/3 1 5/3]);
set(hc,'YTickLabel',{'Exterior','Boundary','Intertior'});
axis equal; view(3); axis tight;  grid on;  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,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,V_hex,[],5);
logicUni=F_count==1; %Logic for boundary faces

Fq=Fp(logicUni,:);
Vq=Vp;
[Fq,Vq,~]=patchCleanUnused(Fq,Vq);

```

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;
patch('Faces',Fs,'Vertices',Vs,'FaceColor','g','EdgeColor','none','FaceAlpha',0.25);
patch('Faces',Ft,'Vertices',Vt,'FaceColor','b','EdgeColor','none','FaceAlpha',0.5);

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);

[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.2AaYQV';

modelName='tetGenModel';
smeshName=[modelName,'.smesh'];

smeshStruct.stringOpt=stringOpt;
smeshStruct.Faces=F;
smeshStruct.Nodes=V;
smeshStruct.holePoints=V_holes;
smeshStruct.faceBoundaryMarker=faceBoundaryMarker; %Face boundary markers
smeshStruct.regionPoints=V_regions; %region points
smeshStruct.regionA=regionA;
smeshStruct.minRegionMarker=2; %Minimum region marker
smeshStruct.smeshName=smeshName;
```

Mesh model using tetrahedral elements using tetGen (see: http://wias-berlin.de/software/tetgen/)

```[meshOutput]=runTetGen(smeshStruct); %Run tetGen
```
```
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- TETGEN Tetrahedral meshing --- 11-Dec-2020 12:16:07
Warning: smeshStruct.smeshName input will be replaced by smeshStruct.modelName
in future releases!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing SMESH file --- 11-Dec-2020 12:16:07
Warning: smeshStruct.smeshName input will be replaced by smeshStruct.modelName
in future releases!
--- Done --- 11-Dec-2020 12:16:07
--- Running TetGen to mesh input boundary--- 11-Dec-2020 12:16:07
Opening /mnt/data/MATLAB/GIBBON/data/temp/tetGenModel.smesh.
Initializing memorypools.
tetrahedron per block: 8188.
Size of a point: 128 bytes.
Size of a tetrahedron: 112 (112) bytes.
Size of a shellface: 192 (192) bytes.
Initializing robust predicates.
sizeof(double) =  8
machine epsilon =   2.22045e-16 [IEEE 754 64-bit macheps]
Permuting vertices.
Sorting vertices.
Incrementally inserting vertices.
7188 (10775) subfaces (segments).
Recovering segments.
10768 (7) segments are recovered (missing).
10773 (2) segments are recovered (missing).
Added 2 Steiner points in volume.
Recovering facets.
7188 (0) subfaces are recovered (missing).
Removed 2 exterior Steiner vertices.
Found 1 domain.
Initial obj = 3287.0271812079686
Final obj  = 3286.9904816687053
Min dihedral   angle  = 0.
Performed 1 flips.
Optimization level  = 2.
Optimization scheme = 7.
Number of iteration = 3.
Min_Max dihed angle = 165.
Removed 361 edges.
0 duplicated vertices are removed.
2 unused vertices are removed.
--- Done --- 11-Dec-2020 12:16:08

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Importing TetGen files --- 11-Dec-2020 12:16:08
--- Done --- 11-Dec-2020 12:16:08
```

Access model element and patch data

```F_tet=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,V);

E_tet=ind2(E_tet);
E_hex=ind2(E_hex);
indBoundary=ind2(indBoundary);

E={E_tet, E_hex};
```

PLOTTING MODEL

```hf1=cFigure;
% subplot(1,2,1);
title('XY cut view of mixed mesh','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

%Selecting half of the model to see interior
X=V(:,3); XE=mean(X(E{1}),2);
L=XE<mean(X);
[Fs,Cs]=element2patch(E{1}(L,:),C_tet(L));

%Selecting half of the model to see interior
X=V(:,3); XE=mean(X(E{2}),2);
L=XE<mean(X);
[Fh,Ch]=element2patch(E{2}(L,:),C_hex(L));

hps=patch('Faces',Fs,'Vertices',V,'FaceColor','g','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
hps2=patch('Faces',Fh,'Vertices',V,'FaceColor','b','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
set(gca,'FontSize',fontSize);
drawnow;

hf1=cFigure;
% subplot(1,2,2);
title('XZ cut view of mixed mesh','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{1}),2);
L=XE>mean(X);
[Fs,Cs]=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);
[Fh,Ch]=element2patch(E{2}(L,:),C_hex(L));

hps=patch('Faces',Fs,'Vertices',V,'FaceColor','g','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
hps2=patch('Faces',Fh,'Vertices',V,'FaceColor','b','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);

view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
set(gca,'FontSize',fontSize);

drawnow;
```
```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);
```

PLOTTING MODEL

```hf1=cFigure;
% subplot(1,2,1);
title('XY cut view of mixed mesh','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

%Selecting half of the model to see interior
X=V(:,3); XE=mean(X(E{1}),2);
L=XE<mean(X);
[Fs,Cs]=element2patch(E{1}(L,:),C_tet(L));

%Selecting half of the model to see interior
X=V(:,3); XE=mean(X(E{2}),2);
L=XE<mean(X);
[Fh,Ch]=element2patch(E{2}(L,:),C_hex(L));

hps=patch('Faces',Fs,'Vertices',VS,'FaceColor','g','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
hps2=patch('Faces',Fh,'Vertices',VS,'FaceColor','b','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
set(gca,'FontSize',fontSize);
drawnow;

hf1=cFigure;
% subplot(1,2,2);
title('XZ cut view of mixed mesh','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

%Selecting half of the model to see interior
X=V(:,2); XE=mean(X(E{1}),2);
L=XE>mean(X);
[Fs,Cs]=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);
[Fh,Ch]=element2patch(E{2}(L,:),C_hex(L));

hps=patch('Faces',Fs,'Vertices',VS,'FaceColor','g','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);
hps2=patch('Faces',Fh,'Vertices',VS,'FaceColor','b','lineWidth',edgeWidth,'edgeColor',edgeColor,'FaceAlpha',1);

view(3); axis tight;  axis equal;  grid on;
colormap(autumn);
set(gca,'FontSize',fontSize);

drawnow;
```

GIBBON

Kevin M. Moerman (kevinmoerman@hotmail.com)

GIBBON footer text

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/.