DEMO_febio_0060_vertebrae_disc_01

Below is a demonstration for:

Contents

Keywords

clear; close all; clc;

Plot settings

fontSize=15;
faceAlpha1=1;
faceAlpha2=0.3;
markerSize1=25;
markerSize2=10;
lineWidth=2;
boneColor=[1.0000    0.9500    0.8000];
discColor=[0.8500    0.5000    0.3000];

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=fullfile(savePath,[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
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_energy_out.txt']; %Log file name for exporting strain energy density

%Geometric parameters
vertebraOffset=[0 0 30];
discHeight=12;
rotAngleVert2=0;
volumeFactorDisc=1;
annulusFibrosusFraction=0.5;
numLayers=3;
% Material parameters
% Disc
c1_1=1; %Shear-modulus-like parameter in MPa
m1_1=2; %Material parameter setting degree of non-linearity
k_1=c1_1*100; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=10; %Optimum number of iterations
max_retries=8; %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-20;
runMode='external';

%Boundary condition parameters
forceApplied=1;
displacementMagnitude=-1;

Get vertebra bone model

[F1,V1]=graphicsModels('vertebra');
rotAngle=-6;
R=euler2DCM([0 (rotAngle./180)*pi 0]);
V1=V1*R;
pointSpacing=mean(patchEdgeLengths(F1,V1));
V1F=patchCentre(F1,V1);

N=patchNormal(F1,V1);
nz=[0 0 1];
a=acos(dot(N,nz(ones(size(N,1),1),:),2));

D=sqrt(sum(V1F(:,[1 2]).^2,2));
logicTop=D<24 & a<0.27;
indNotInLogic=unique(F1(~logicTop,:));
logicTop=logicTop & ~any(ismember(F1,indNotInLogic),2);
logicTop=triSurfLogicSharpFix(F1,logicTop,3);

D=minDist(V1F(:,[1 2]),V1F(logicTop,[1 2]));
logicBottom=D<3 & a>pi-0.4;
logicBottom=triSurfLogicSharpFix(F1,logicBottom,3);

indNotInLogic=unique(F1(~logicBottom,:));
logicBottom=logicBottom & ~any(ismember(F1,indNotInLogic),2);
logicBottom=triSurfLogicSharpFix(F1,logicBottom,3);

C1=zeros(size(F1,1),1);
C1(logicTop)=1;
C1(logicBottom)=2;
Eb1=patchBoundary(F1(C1==1,:),V1);
indList1=edgeListToCurve(Eb1);
indList1=indList1(1:end-1);

Eb2=patchBoundary(F1(C1==2,:),V1);
indList2=edgeListToCurve(Eb2);
indList2=indList2(1:end-1);
n=[numel(indList1) numel(indList2)];
[nMax,indMax]=max(n);
nn=abs(diff(n));
switch indMax
    case 1
        [F1,V1,Eb2,C1]=triSurfSplitBoundary(F1,V1,Eb2,size(Eb2,1)+nn,C1);
        indList2=edgeListToCurve(Eb2);
        indList2=indList2(1:end-1);
    case 2
        [F1,V1,Eb1,C1]=triSurfSplitBoundary(F1,V1,Eb1,size(Eb1,1)+nn,C1);
        indList1=edgeListToCurve(Eb1);
        indList1=indList1(1:end-1);
end

indList2=fliplr(indList2);
nc=10;
ns=50;
[V1]=smoothCurve(F1,V1,nc,ns,indList1);
[V1]=smoothCurve(F1,V1,nc,ns,indList2);

P11=mean(V1(indList1,:),1);
N11=mean(patchNormal(F1(C1==1,:),V1),1);

P12=mean(V1(indList2,:),1);
N12=-mean(patchNormal(F1(C1==2,:),V1),1);
cFigure; hold on;
gpatch(F1,V1,C1,'k',1);
plotV(V1(indList1,:),'r.-','MarkerSize',25,'LineWidth',3);
plotV(V1(indList2,:),'r.-','MarkerSize',25,'LineWidth',3);
% plotV(V1(indRigid,:),'c.','MarkerSize',15);
quiverVec(P11,N11,5,'r');
quiverVec(P12,N12,5,'r');

axisGeom(gca,fontSize);
colormap gjet; icolorbar;
camlight('headlight');
gdrawnow;
F2=F1;
V2=V1;
C2=C1;

R2=euler2DCM([0 (rotAngleVert2./180)*pi 0]);
V2=V2*R2;
V2=V2+vertebraOffset(ones(size(V1,1),1),:);
V2(:,3)=V2(:,3)+discHeight;

P21=mean(V2(indList1,:),1);
N21=mean(patchNormal(F2(C2==1,:),V2),1);

P22=mean(V2(indList2,:),1);
N22=-mean(patchNormal(F2(C2==2,:),V2),1);
f=discHeight/3;
p=[P11;P11+f*N11; P22-f*N22; P22];

numStepsCurve=ceil(discHeight./pointSpacing);
Vg=bezierCurve(p,numStepsCurve);

[Fds,Vds,Cds]=sweepLoft(V1(indList1,:),V2(indList2,:),N11,N12,Vg,numStepsCurve);
[~,~,Nd]=patchNormal(Fds,Vds);
CVds=faceToVertexMeasure(Fds,Vds,Cds);
CVds=CVds-min(CVds(:));
CVds=CVds./max(CVds(:));
CVds=abs(CVds-0.5);
CVds=CVds./max(CVds(:));
CVds=CVds.^2;
CVds=1-CVds;
Vds=Vds+discHeight/10*Nd.*CVds;

[Fds,Vds]=quad2tri(Fds,Vds,'a');

clear cPar;
cPar.n=5;
cPar.Method='HC';
cPar.RigidConstraints=unique(patchBoundary(Fds,Vds));
[Vds]=patchSmooth(Fds,Vds,[],cPar);

[Fd,Vd,Cd]=joinElementSets({Fds,fliplr(F2(C2==2,:)),fliplr(F1(C1==1,:))},{Vds,V2,V1});
[Fd,Vd]=patchCleanUnused(Fd,Vd);
[Fd,Vd]=mergeVertices(Fd,Vd);

Visualize imported surfaces

cFigure; hold on;
gpatch(F1,V1,'rw','none',0.5);
gpatch(F2,V2,'gw','none',0.5);
gpatch(Fd,Vd,Cd,'k',1);

% plotV(Vg,'r.-','MarkerSize',25,'LineWidth',3);

camlight('headlight');
axisGeom(gca,fontSize);
colormap gjet; icolorbar;
gdrawnow;

Mesh disc with tetrahedral elements

Tet meshing is based on tetgen. TetGen requires a interior points for regions to be meshed, as well as intertior points for holes.

Define region points

[V_region]=getInnerPoint(Fd,Vd);

Visualize interior points

cFigure; hold on;
hp1=gpatch(Fd,Vd,'kw','none',0.5);
patchNormPlot(Fd,Vd);
hp2=plotV(V_region,'r.','markerSize',markerSize1);
legend([hp1 hp2],{'Disc mesh','Region point'});
axisGeom(gca,fontSize);
camlight('headlight');
gdrawnow;

Mesh using tetgen

inputStruct.stringOpt='-pq1.2AaY'; %TetGen option string
inputStruct.Faces=Fd; %The faces
inputStruct.Nodes=Vd; %The vertices
inputStruct.holePoints=[]; %The hole interior points
inputStruct.faceBoundaryMarker=Cd; %Face boundary markers
inputStruct.regionPoints=V_region; %The region interior points
inputStruct.regionA=tetVolMeanEst(Fd,Vd)*volumeFactorDisc; %Volume for regular tets

Mesh model using tetrahedral elements using tetGen

[meshOutput]=runTetGen(inputStruct); %Run tetGen
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- TETGEN Tetrahedral meshing --- 07-Jul-2020 10:03:14
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing SMESH file --- 07-Jul-2020 10:03:14
----> Adding node field
----> Adding facet field
----> Adding holes specification
----> Adding region specification
--- Done --- 07-Jul-2020 10:03:14
--- Running TetGen to mesh input boundary--- 07-Jul-2020 10:03:14
Opening /mnt/data/MATLAB/GIBBON/data/temp/temp.smesh.
Delaunizing vertices...
Delaunay seconds:  0.018735
Creating surface mesh ...
Surface mesh seconds:  0.002958
Recovering boundaries...
Boundary recovery seconds:  0.004153
Removing exterior tetrahedra ...
Spreading region attributes.
Exterior tets removal seconds:  0.002478
Recovering Delaunayness...
Delaunay recovery seconds:  0.002329
Refining mesh...
Refinement seconds:  0.169956
Optimizing mesh...
Optimization seconds:  0.013369

Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.node.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.ele.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.face.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.edge.

Output seconds:  0.156609
Total running seconds:  0.370885

Statistics:

  Input points: 2066
  Input facets: 4128
  Input segments: 6192
  Input holes: 0
  Input regions: 1

  Mesh points: 12558
  Mesh tetrahedra: 73607
  Mesh faces: 149278
  Mesh faces on exterior boundary: 4128
  Mesh faces on input facets: 4128
  Mesh edges on input segments: 6192
  Steiner points inside domain: 10492

--- Done --- 07-Jul-2020 10:03:14
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Importing TetGen files --- 07-Jul-2020 10:03:14
--- Done --- 07-Jul-2020 10:03:15

Access model element and patch data

Fb_disc=meshOutput.facesBoundary; %Boundary faces of the disc
Cb_disc=meshOutput.boundaryMarker; %Boundary marker/color data for the disc
V_disc=meshOutput.nodes; %The vertices/nodes
E_disc=meshOutput.elements; %The tet4 elements

Visualizing mesh using meshView, see also anim8

meshView(meshOutput);
VE_disc=patchCentre(E_disc,V_disc);
Eb=patchBoundary(Fb_disc(ismember(Cb_disc,[2,3]),:),V_disc);
indB=unique(Fb_disc(ismember(Cb_disc,1),:));
[D,~]=minDist(VE_disc(:,[1 2]),V_disc(indB,[1 2]));
D=D./max(D(:));
D=1-D;
logicAnnulusFibrosus=D<=annulusFibrosusFraction;
D(logicAnnulusFibrosus)=annulusFibrosusFraction;
D=D-min(D(:));
D=D./max(D(:));


[FE,CF]=element2patch(E_disc,logicAnnulusFibrosus);
cFigure;  hold on;
gpatch(FE,V_disc,CF,'k',1);

axisGeom;
camlight headlight;
view(2);
colormap gjet; colorbar;
gdrawnow;
VEb=patchCentre(Eb,V_disc);
t = atan2(VEb(:,2),VEb(:,1));
T = atan2(VE_disc(:,2),VE_disc(:,1));


VEb=patchCentre(Eb,V_disc);
NE=vecnormalize(V_disc(Eb(:,1),:)-V_disc(Eb(:,2),:));
NE(:,3)=0;
NE=vecnormalize(NE);
[~,indMin]=minDist(T,t);
% z=[0 0 1];
% NFC=cross(NF,z(ones(size(NF,1),1),:),2);
NE_disc=NE(indMin,:);

% NE_disc(:,[1 2])=NE_disc(:,[1 2]).*sin(numLayers.*D(:,ones(1,2)).*2*pi);
% NE_disc(:,3)=cos(numLayers.*D.*2*pi);
% NE_disc(logicAnnulusFibrosus,:)=NaN;

logicEven=iseven(round(D.*numLayers));
NE_disc(logicEven,3)=1;
NE_disc(~logicEven,3)=-1;
NE_disc=vecnormalize(NE_disc);
NE_disc(logicAnnulusFibrosus,:)=NaN;

cFigure;  hold on;
gpatch(Fb_disc,V_disc,'w','none',0.1);
gpatch(Eb,V_disc,'none','k',1,3);
% plotV(VE_disc,'k.','MarkerSize',1);
plotV(mean(V_disc,1),'r.','MarkerSize',25);
quiverVec(VE_disc(1:10:end,:),NE_disc(1:10:end,:),1,D(1:10:end,:));
axisGeom;
camlight headlight;
view(2);
colormap gjet;
gdrawnow;

Joining node sets

V=[V_disc;V1;V2]; %Combined node sets
E1=F1+size(V_disc,1); %Fixed element indices
E2=F2+size(V_disc,1)+size(V1,1); %Fixed element indices

numDigitsMerge=6-numOrder(pointSpacing);
[~,indKeep,indFix]=unique(pround(V,numDigitsMerge),'rows');
V=V(indKeep,:);
E_disc=indFix(E_disc);
E1=indFix(E1);
E2=indFix(E2);
Fb_disc=indFix(Fb_disc);

Define boundary conditions

Visualize BC's

cFigure; hold on;
title('Boundary conditions');
hp(1)=gpatch(Fb_disc,V,'w','none',0.5);
hp(2)=gpatch(E1,V,'rw','none',1);
hp(3)=gpatch(E2,V,'bw','none',1);
legend(hp,{'Disc','Constrained vertebra','BC vertebra'});

axisGeom;
camlight headlight;
gdrawnow;

Defining the FEBio input structure

See also febioStructTemplate and febioStruct2xml and the FEBio user manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='2.5';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis.ATTR.type='static';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;
febio_spec.Control.max_refs=max_refs;
febio_spec.Control.max_ups=max_ups;
febio_spec.Control.symmetric_stiffness=symmetric_stiffness;
febio_spec.Control.min_residual=min_residual;

%Material section
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1_1;
febio_spec.Material.material{1}.m1=m1_1;
febio_spec.Material.material{1}.c2=c1_1;
febio_spec.Material.material{1}.m2=-m1_1;
febio_spec.Material.material{1}.k=k_1;

febio_spec.Material.material{2}.ATTR.type='rigid body';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.density=1;
febio_spec.Material.material{2}.center_of_mass=mean(V1);

febio_spec.Material.material{3}.ATTR.type='rigid body';
febio_spec.Material.material{3}.ATTR.id=3;
febio_spec.Material.material{3}.density=1;
febio_spec.Material.material{3}.center_of_mass=mean(V2);

%Geometry section
% -> Nodes
febio_spec.Geometry.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Geometry.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Geometry.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
febio_spec.Geometry.Elements{1}.ATTR.type='tet4'; %Element type of this set
febio_spec.Geometry.Elements{1}.ATTR.mat=1; %material index for this set
febio_spec.Geometry.Elements{1}.ATTR.name='disc'; %Name of the element set
febio_spec.Geometry.Elements{1}.elem.ATTR.id=(1:1:size(E_disc,1))'; %Element id's
febio_spec.Geometry.Elements{1}.elem.VAL=E_disc;

febio_spec.Geometry.Elements{2}.ATTR.type='tri3'; %Element type of this set
febio_spec.Geometry.Elements{2}.ATTR.mat=2; %material index for this set
febio_spec.Geometry.Elements{2}.ATTR.name='Bone1'; %Name of the element set
febio_spec.Geometry.Elements{2}.elem.ATTR.id=size(E_disc,1)+(1:1:size(E1,1))'; %Element id's
febio_spec.Geometry.Elements{2}.elem.VAL=E1;

febio_spec.Geometry.Elements{3}.ATTR.type='tri3'; %Element type of this set
febio_spec.Geometry.Elements{3}.ATTR.mat=3; %material index for this set
febio_spec.Geometry.Elements{3}.ATTR.name='Bone2'; %Name of the element set
febio_spec.Geometry.Elements{3}.elem.ATTR.id=size(E_disc,1)+size(E1,1)+(1:1:size(E2,1))'; %Element id's
febio_spec.Geometry.Elements{3}.elem.VAL=E2;

% -> NodeSets
febio_spec.Geometry.NodeSet{1}.ATTR.name='All';
febio_spec.Geometry.NodeSet{1}.node.ATTR.id=(1:1:size(V,1))';

%Boundary condition section
% -> Prescribed boundary conditions on the rigid body
febio_spec.Boundary.rigid_body{1}.ATTR.mat=2;
febio_spec.Boundary.rigid_body{1}.fixed{1}.ATTR.bc='x';
febio_spec.Boundary.rigid_body{1}.fixed{2}.ATTR.bc='y';
febio_spec.Boundary.rigid_body{1}.fixed{3}.ATTR.bc='z';
febio_spec.Boundary.rigid_body{1}.fixed{4}.ATTR.bc='Rx';
febio_spec.Boundary.rigid_body{1}.fixed{5}.ATTR.bc='Ry';
febio_spec.Boundary.rigid_body{1}.fixed{6}.ATTR.bc='Rz';

febio_spec.Boundary.rigid_body{2}.ATTR.mat=3;
febio_spec.Boundary.rigid_body{2}.fixed{1}.ATTR.bc='x';
febio_spec.Boundary.rigid_body{2}.fixed{2}.ATTR.bc='y';
febio_spec.Boundary.rigid_body{2}.fixed{3}.ATTR.bc='Rx';
febio_spec.Boundary.rigid_body{2}.fixed{4}.ATTR.bc='Ry';
febio_spec.Boundary.rigid_body{2}.fixed{5}.ATTR.bc='Rz';

febio_spec.Boundary.rigid_body{2}.prescribed.ATTR.bc='z';
febio_spec.Boundary.rigid_body{2}.prescribed.ATTR.lc=1;
febio_spec.Boundary.rigid_body{2}.prescribed.VAL=displacementMagnitude;

% febio_spec.Boundary.rigid_body{2}.force.ATTR.bc='z';
% febio_spec.Boundary.rigid_body{2}.force.ATTR.lc=1;
% febio_spec.Boundary.rigid_body{2}.force.VAL=forceBody;

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.rigid_body_data{1}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.rigid_body_data{1}.ATTR.data='Fx;Fy;Fz';
febio_spec.Output.logfile.rigid_body_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.rigid_body_data{1}.VAL=2; %Rigid body material id

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E_disc,1);

Running the FEBio analysis

To run the analysis defined by the created FEBio input file the runMonitorFEBio function is used. The input for this function is a structure defining job settings e.g. the FEBio input file name. The optional output runFlag informs the user if the analysis was run succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode=runMode;%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

Export input file

Exporting the febio_spec structure to an FEBio input file is done using the febioStruct2xml function.

    febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

Run febio

    [runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- STARTING FEBIO JOB --- 07-Jul-2020 10:03:26
Waiting for log file...
Proceeding to check log file...07-Jul-2020 10:03:27
------- converged at time : 0.1
------- converged at time : 0.2
------- converged at time : 0.3
------- converged at time : 0.4
------- converged at time : 0.5
------- converged at time : 0.6
------- converged at time : 0.7
------- converged at time : 0.8
------- converged at time : 0.9
------- converged at time : 1
--- Done --- 07-Jul-2020 10:04:16

Importing rigid body reaction forces from a log file

        [time_mat, R_force_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_force)); %Nodal forces
        time_mat=[0; time_mat(:)]; %Time

        F_reaction=[0 0 0; R_force_mat(:,2:end)]; %Get force, add zeros, remove rigidbody id column

        Fz_final=F_reaction(end,3);

Importing nodal displacements from a log file

[time_mat, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements
time_mat=[0; time_mat(:)]; %Time

N_disp_mat=N_disp_mat(:,2:end,:);
sizImport=size(N_disp_mat);
sizImport(3)=sizImport(3)+1;
N_disp_mat_n=zeros(sizImport);
N_disp_mat_n(:,:,2:end)=N_disp_mat;
N_disp_mat=N_disp_mat_n;
DN=N_disp_mat(:,:,end);

V_def=V+DN;
V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
X_DEF=V_DEF(:,1,:);
Y_DEF=V_DEF(:,2,:);
Z_DEF=V_DEF(:,3,:);

C=sqrt(sum(DN(:,3).^2,2));
[CF]=vertexToFaceMeasure(Fb_disc,C);

Importing element strain energies from a log file

[~,E_energy,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy)); %Element strain energy

%Remove nodal index column
E_energy=E_energy(:,2:end,:);

%Add initial state i.e. zero energy
sizImport=size(E_energy);
sizImport(3)=sizImport(3)+1;
E_energy_mat_n=zeros(sizImport);
E_energy_mat_n(:,:,2:end)=E_energy;
E_energy=E_energy_mat_n;

[FE_disc,C_energy_disc]=element2patch(E_disc,E_energy(:,:,end));
%     [FE_disc,C_energy_disc]=element2patch(E_disc,E_energy(1:size(E_disc,1),:,1));
indBoundaryFacesdisc=tesBoundary(FE_disc,V);

Plot animation

Plotting the simulated results using anim8 to visualize and animate deformations

% Create basic view and store graphics handle to initiate animation
hf=cFigure; hold on; %Open figure
gtitle([febioFebFileNamePart,': Press play to animate']);
CV=faceToVertexMeasure(FE_disc(indBoundaryFacesdisc,:),V_def,C_energy_disc(indBoundaryFacesdisc));
hp1=gpatch(Fb_disc,V_def,CV,'none',1); %Add graphics object to animate
hp1.FaceColor='Interp';
hp2=gpatch([E1;E2],V_def,'w','none',0.5); %Add graphics object to animate

axisGeom(gca,fontSize);
colormap(gjet(250)); colorbar;
caxis([0 max(C_energy_disc)/25]);
axis([min(X_DEF(:)) max(X_DEF(:)) min(Y_DEF(:)) max(Y_DEF(:)) min(Z_DEF(:)) max(Z_DEF(:))]);
camlight headlight;

% Set up animation features
animStruct.Time=time_mat; %The time vector
for qt=1:1:size(N_disp_mat,3) %Loop over time increments
    DN=N_disp_mat(:,:,qt); %Current displacement
    V_def=V+DN; %Current nodal coordinates

    [FE_disc,C_energy_disc]=element2patch(E_disc,E_energy(:,:,qt));
    CV=faceToVertexMeasure(FE_disc(indBoundaryFacesdisc,:),V_def,C_energy_disc(indBoundaryFacesdisc));

    %Set entries in animation structure
    animStruct.Handles{qt}=[hp1 hp1 hp2]; %Handles of objects to animate
    animStruct.Props{qt}={'Vertices','CData','Vertices'}; %Properties of objects to animate
    animStruct.Set{qt}={V_def,CV,V_def}; %Property values for to set in order to animate
end
anim8(hf,animStruct); %Initiate animation feature
gdrawnow;
function [V1]=smoothCurve(F1,V1,nc,ns,indList1)
clear cPar
cPar.n=nc;
cPar.Method='HC';
e=[(1:numel(indList1))' [(2:numel(indList1))';1]];
v=V1(indList1,:);
[v]=patchSmooth(e,v,[],cPar);
V1(indList1,:)=v;

indTouch=unique(F1(any(ismember(F1,indList1),2),:));
indTouch=unique(F1(any(ismember(F1,indTouch),2),:));
logicRigid=true(size(V1,1),1);
logicRigid(indTouch)=0;
logicRigid(indList1)=1;
indRigid=find(logicRigid);

cPar.n=ns;
cPar.Method='HC';
cPar.RigidConstraints=indRigid;
[V1]=patchSmooth(F1,V1,[],cPar);
end

GIBBON www.gibboncode.org

Kevin Mattheus Moerman, [email protected]

GIBBON discer 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) 2019 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/.

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