You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: DEMO_febio_0044_mammography_01.m
    path of file: ./repos/GIBBON/docs
    
    the code of the file until where you have to start completion: %% DEMO_febio_0044_mammography_01.m
% Below is a demonstration for:
%
% * Building geometry for a slab with hexahedral elements, and a
% triangulated sphere.
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 4.0
% * febio, FEBio
% * indentation
% * contact, sliding, sticky, friction
% * rigid body constraints
% * tetrahedral elements, tet4
% * triangular elements, tri3
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
markerSize1=40;
markerSize2=25;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force

%Breast geometry
r_breastNominal=40; %Breast radius
sy=1; %Y-direction radial scaling to force deviation from hemi-sphere
r1=r_breastNominal/2.5; %Radius alveola
r2=r_breastNominal/7; %Radius nipple
rm=mean([r1 r2]); %radius between outer alveola radius and outer nipple radius
w=(r1-r2)/20; %Height of alveola
h=r2; %Height of alveola
gravityShiftPercentage=0.5; %Percentage shift due to gravity
dx=r_breastNominal.*gravityShiftPercentage; %X-shift due to gravity

%Loading
quasiCompressionLevel=0.3; %Sets displacement as function of nominal breast radius
loadAngle=(0/180)*pi;
distanceChestPlate=12; %r_breastNominal/7

%Meshing parameters
nRefine=3; %Number of refine steps for hemi-sphere
volumeFactor=3;
plateHeight=r_breastNominal;

%Material parameter set
k_factor=50; %Bulk modulus factor

c1_1=1e-3; %Shear-modulus-like parameter
m1_1=2; %Material parameter setting degree of non-linearity
k_1=c1_1*k_factor; %Bulk modulus

c1_2=c1_1*5; %Shear-modulus-like parameter
m1_2=m1_1; %Material parameter setting degree of non-linearity
k_2=c1_2*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
symmetric_stiffness=0;
min_residual=1e-18;

%Set run mode
runMode='external'; %'internal'

%Contact parameters
contactPenalty=15;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.1;

%% Create hemi-sphere
[F,V,C_hemiSphereLabel]=hemiSphereMesh(nRefine,r_breastNominal,1); %Construct hemi-shere mesh

%% Change shape of hemi-sphere to create basic breast model

%scale Y direction to deviate from hemi-sphere. 
V(:,2)=V(:,2)*sy;

indExclude=unique(F(C_hemiSphereLabel==2,:));
logicExclude=false(size(V,1),1);
logicExclude(indExclude)=1;

dt=sqrt(sum(V(:,[1 2]).^2,2));
logicHigh1=dt<r1 & ~logicExclude;
logicHigh2=dt<r2 & ~logicExclude;

C_skin=double(logicHigh1);
C_skin(logicHigh2)=2;

t=linspace(0,2*pi,500);
x=rm*sin(t);
y=rm*cos(t);
vc=[x(:) y(:)];

[d]=minDist(V(C_skin==1,[1 2]),vc);
dtt=d.^3;
dtt=dtt-min(dtt);
dtt=dtt./max(dtt);
dtt=abs(dtt-1)*w;
V(C_skin==1,3)=V(C_skin==1,3)+dtt;

f=V(:,3); %Z-coordinates [min(Z) max(Z)]
f=f-min(f(:)); %Z-coordinates [0 max(Z)-min(Z)]
f=f./max(f(:)); % [0 1] %Normalized Z-coordinate parameter
V(:,1)=V(:,1)+dx.*f;

dtt=dt(C_skin==2).^3;
dtt=dtt-min(dtt);
dtt=dtt./max(dtt);
dtt=abs(dtt-1)*h;
V(C_skin==2,3)=V(C_skin==2,3)+dtt;

%%

cFigure; hold on;
gpatch(F,V,'w','none',1);
axisGeom;
camlight headlight;
gdrawnow;

%%

breastVolume=triSurfVolume(F,V)';

%% Rotate model
R=euler2DCM([pi -0.5*pi 0
    
function of nominal breast radius
loadAngle=(0/180)*pi;
distanceChestPlate=12; %r_breastNominal/7

%Meshing parameters
nRefine=3; %Number of refine steps for hemi-sphere
volumeFactor=3;
plateHeight=r_breastNominal;

%Material parameter set
k_factor=50; %Bulk modulus factor

c1_1=1e-3; %Shear-modulus-like parameter
m1_1=2; %Material parameter setting degree of non-linearity
k_1=c1_1*k_factor; %Bulk modulus

c1_2=c1_1*5; %Shear-modulus-like parameter
m1_2=m1_1; %Material parameter setting degree of non-linearity
k_2=c1_2*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
symmetric_stiffness=0;
min_residual=1e-18;

%Set run mode
runMode='external'; %'internal'

%Contact parameters
contactPenalty=15;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.1;

%% Create hemi-sphere
[F,V,C_hemiSphereLabel]=hemiSphereMesh(nRefine,r_breastNominal,1); %Construct hemi-shere mesh

%% Change shape of hemi-sphere to create basic breast model

%scale Y direction to deviate from hemi-sphere. 
V(:,2)=V(:,2)*sy;

indExclude=unique(F(C_hemiSphereLabel==2,:));
logicExclude=false(size(V,1),1);
logicExclude(indExclude)=1;

dt=sqrt(sum(V(:,[1 2]).^2,2));
logicHigh1=dt<r1 & ~logicExclude;
logicHigh2=dt<r2 & ~logicExclude;

C_skin=double(logicHigh1);
C_skin(logicHigh2)=2;

t=linspace(0,2*pi,500);
x=rm*sin(t);
y=rm*cos(t);
vc=[x(:) y(:)];

[d]=minDist(V(C_skin==1,[1 2]),vc);
dtt=d.^3;
dtt=dtt-min(dtt);
dtt=dtt./max(dtt);
dtt=abs(dtt-1)*w;
V(C_skin==1,3)=V(C_skin==1,3)+dtt;

f=V(:,3); %Z-coordinates [min(Z) max(Z)]
f=f-min(f(:)); %Z-coordinates [0 max(Z)-min(Z)]
f=f./max(f(:)); % [0 1] %Normalized Z-coordinate parameter
V(:,1)=V(:,1)+dx.*f;

dtt=dt(C_skin==2).^3;
dtt=dtt-min(dtt);
dtt=dtt./max(dtt);
dtt=abs(dtt-1)*h;
V(C_skin==2,3)=V(C_skin==2,3)+dtt;

%%

cFigure; hold on;
gpatch(F,V,'w','none',1);
axisGeom;
camlight headlight;
gdrawnow;

%%

breastVolume=triSurfVolume(F,V)';

%% Rotate model
R=euler2DCM([pi -0.5*pi 0]);
V=V*R;

%%

[Fs,Vs]=geoSphere(2,r_breastNominal/6);
Vs(:,1)=Vs(:,1)-r_breastNominal/2;
Vs(:,2)=Vs(:,2)-r_breastNominal/4;
Vs(:,3)=Vs(:,3)-r_breastNominal/2;

%%

% cFigure; hold on;
% gpatch(F,V,C_skin,'none',1);
% axisGeom;
% camlight headlight;
% icolorbar;
% gdrawnow;

%%

C=[C_hemiSphereLabel;(max(C_hemiSphereLabel(:))+1)*ones(size(Fs,1),1)];
F=[F;Fs+size(V,1)];
V=[V;Vs];

%%

cFigure; hold on;
gpatch(F,V,C,'none',1);
axisGeom;
camlight headlight;
icolorbar;
gdrawnow;

%% Get interior points
[V_in1]=getInnerPoint({F,Fs},{V,Vs});
[V_in2]=getInnerPoint(Fs,Vs);

V_regions=[V_in1; V_in2];

cFigure; hold on;
gpatch(F,V,C,'none',0.5);
% gpatch(Fs,Vs,'kw','none',1);
plotV(V_regions,'k.','markerSize',50);
axisGeom;
camlight headlight;
gdrawnow;

%%

faceBoundaryMarker=C;

[regionA]=tetVolMeanEst(F,V); %Volume for regular tets
inputStruct.stringOpt='-pq1.2AaY';
inputStruct.Faces=fliplr(F);
inputStruct.Nodes=V;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=faceBoundaryMarker; %Face boundary markers
inputStruct.regionPoints=V_regions; %region points
inputStruct.regionA=regionA*ones(size(V_regions,1),1)*volumeFactor;
inputStruct.minRegionMarker=2; %Minimum region marker

% Mesh model using tetrahedral elements using tetGen
[meshOutput]=runTetGen(inputStruct); %Run tetGen

% Access model element and patch data
Fb=meshOutput.facesBoundary;
Cb=meshOutput.boundaryMarker;
V=meshOutput.nodes;
CE=meshOutput.elementMaterialID;
E=meshOutput.elements;

%%

Q=euler2DCM([loadAngle 0 0]);
V=V*Q; 

%%

cFigure; hold on;
hp=gpatch(Fb,V,Cb,'none',0.5);
axisGeom;
camlight headlight;
icolorbar;
gdrawnow;

%% Visualizing mesh using |meshView|, see also |anim8|

meshView(meshOutput);

%% Split element sets

E1=E(meshOutput.elementMaterialID==-2,:);
E2=E(meshOutput.elementMaterialID==-3,:);

%% Building plate models

%Basic side curve of plate
Vt=[0 0 0; 0 0 plateHeight; 1.5*r_breastNominal 0 plateHeight; 1.5*r_breastNominal 0 0; ];

%Fillet side curve of plate
rFillet=6; %Fillet radius
np=25; %Number of points used to construct each fillet edge
closedLoopOption=0; %Use 1 if curve represents a closed loop but containes unique points
[Vc]=filletCurve(Vt,rFillet,np,closedLoopOption);

pointSpacingPlate=mean(patchEdgeLengths(Fb,V))/2; % Get point spacing from mesh

%Extrude to form plate
cPar.pointSpacing=pointSpacingPlate;
cPar.depth=3*(r_breastNominal*sy);
cPar.patchType='quad';
cPar.dir=0;
cPar.closeLoopOpt=0;
[Fp1,Vp1]=polyExtrude(Vc,cPar);
Np1=mean(patchNormal(Fp1,Vp1),1);
if Np1(3)<0
    Fp1=fliplr(Fp1);
end