DEMO_febio_0054_lattice_hydrostatic_01

Below is a demonstration for:

Contents

Keywords

clear; close all; clc;

Plot settings

Plot settings

fontSize=15;
faceAlpha1=0.8;
faceAlpha2=1;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
markerSize=25;
cMap=gjet(4);

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_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress
febioLogFileName_stiffness=[febioFebFileNamePart,'_stiffness_out.txt']; %Log file name for exporting stiffness

%Specifying dimensions and number of elements
sampleSize=10;
latticeType=1;

%Define applied displacement
J_final=0.7; %Final Jacobian or volume ration
lambdaFinal=J_final^(1/3); %Stretch values in all directions
displacementMagnitude=((lambdaFinal*sampleSize)-sampleSize)/2; %The displacement magnitude

%Material parameter set
c1=1; %Shear-modulus-like parameter
m1=2;
k=50*c1;

% FEA control settings
numTimeSteps=50; %Number of time steps desired
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=25; %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
min_residual=1e-20;
symmetric_stiffness=0;
runMode='external';
%Specifying dimensions and number of elements
r=0.5; %Radii, results in a width of 1
n=4;
nCopies=n*ones(1,3); %Number of offset copies
d=2*r; %Diameter
w=(n-1)*d; %sampleSize

Create lattice

switch latticeType
    case 1 %Octet truss
        [Er,Vr,Cr,Fr,CFr]=rhombicDodecahedronMesh(r,nCopies);
        Vr=Vr./(n-1);
        Vr=Vr*sampleSize;

        [indBoundary]=tesBoundary(Fr,Vr);
        cPar.shrinkFactor=0.15; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=2;
        cPar.latticeSide=2; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [E,V,C]=element2lattice(Er,Vr,cPar); %Get lattice structure

        logicKeep1=~(V(:,1)<=-1e-3);
        logicKeep2=~(V(:,2)<=-1e-3);
        logicKeep3=~(V(:,3)<=-1e-3);
        logicKeep4=~(V(:,1)>=sampleSize+1e-3);
        logicKeep5=~(V(:,2)>=sampleSize+1e-3);
        logicKeep6=~(V(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(E),2)>=4 &...
            sum(logicKeep2(E),2)>=4 &...
            sum(logicKeep3(E),2)>=4 &...
            sum(logicKeep4(E),2)>=4 &...
            sum(logicKeep5(E),2)>=4 &...
            sum(logicKeep6(E),2)>=4;

        E=E(logicKeepEs,:);
        C=C(logicKeepEs,:);
        [E,V,indFix]=patchCleanUnused(E,V);

        % [Es,Vs,~,~]=subHex(Es,Vs,1,1);
        % Cs=repmat(Cs,8,1);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
    case 2 %Rhombic dodecahedron mesh ("dual" of octet truss lattice)
        [Er,Vr,Cr,Fr,CFr]=rhombicDodecahedronMesh(r,nCopies);
        Vr=Vr./(n-1);
        Vr=Vr*sampleSize;

        [indBoundary]=tesBoundary(Fr,Vr);
        cPar.shrinkFactor=0.15; %Strut sides are formed by shrinking the input mesh faces by this factor
        cPar.meshType='hex'; %desired output mesh type
        cPar.indBoundary=indBoundary; %indices of the boundary faces
        cPar.hexSplit=2;
        cPar.latticeSide=1; %1=side 1 the edge lattice, 2=side 2 the dual lattice to the edge lattice
        [E,V,C]=element2lattice(Er,Vr,cPar); %Get lattice structure

        logicKeep1=~(V(:,1)<=-1e-3);
        logicKeep2=~(V(:,2)<=-1e-3);
        logicKeep3=~(V(:,3)<=-1e-3);
        logicKeep4=~(V(:,1)>=sampleSize+1e-3);
        logicKeep5=~(V(:,2)>=sampleSize+1e-3);
        logicKeep6=~(V(:,3)>=sampleSize+1e-3);

        logicKeepEs=sum(logicKeep1(E),2)>=4 &...
            sum(logicKeep2(E),2)>=4 &...
            sum(logicKeep3(E),2)>=4 &...
            sum(logicKeep4(E),2)>=4 &...
            sum(logicKeep5(E),2)>=4 &...
            sum(logicKeep6(E),2)>=4;

        E=E(logicKeepEs,:);
        C=C(logicKeepEs,:);
        [E,V,indFix]=patchCleanUnused(E,V);

        % [Es,Vs,~,~]=subHex(Es,Vs,1,1);
        % Cs=repmat(Cs,8,1);

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        %Get new boundary set
        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
    case 3
        boxDim=sampleSize*[1 1 1];
        boxEl=[2 2 2];
        [meshStruct]=hexMeshBox(boxDim,boxEl);
        Er=meshStruct.E;
        Vr=meshStruct.V;
        minV=min(Vr,[],1); %Get lower left front corner
        Vr=Vr-minV(ones(size(Vr,1),1),:); %Set corner as origin
        [Er,Vr,~]=hex2tet(Er,Vr,[],1); %Convert to tetrahedral elements
        [Fr,Cr]=element2patch(Er,[]); %Patch data for plotting
        [indBoundary]=tesBoundary(Fr,Vr);

        % Create lattice structure
        controlParameter.latticeSide=1;
        controlParameter.numDigitKeep=5; %used for merging nodes
        controlParameter.indBoundary=indBoundary; %indices of the boundary faces
        controlParameter.shrinkFactor=0.15;
        controlParameter.meshType='hex';
        controlParameter.hexSplit=2;

        [E,V,C]=element2lattice(Er,Vr,controlParameter); %Get lattice structure

        % Create patch Data for visualization
        [Fs,CsF]=element2patch(E,C); %Patch data for plotting

        indB=tesBoundary(Fs,V);
        Fb=Fs(indB,:);
end

Visualizing input mesh and lattic structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(Fr,Vr,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,V,'bw','k',1);
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

DEFINE BC's

% Define node set logics
indAll=(1:1:size(V,1))';
logicBoundary=ismember(indAll,Fb);

Z=V(:,3);
logicTop=Z>=(sampleSize-eps(sampleSize))& logicBoundary;
logicBottom=Z<=eps(sampleSize) & logicBoundary;

X=V(:,1);
logicSide1=X>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide2=X<=eps(sampleSize)& logicBoundary;

Y=V(:,2);
logicSide3=Y>=(sampleSize-eps(sampleSize))& logicBoundary;
logicSide4=Y<=eps(sampleSize)& logicBoundary;

%Prescribed force nodes
bcPrescribeListCell{1}=find(logicSide1)';
bcPrescribeListCell{2}=find(logicSide2)';
bcPrescribeListCell{3}=find(logicSide3)';
bcPrescribeListCell{4}=find(logicSide4)';
bcPrescribeListCell{5}=find(logicTop)';
bcPrescribeListCell{6}=find(logicBottom)';

Smoothing lattice

% indKeep=unique([bcPrescribeListCell{:}]);
% [Fb_clean,Vb_clean,indFix]=patchCleanUnused(Fb,Vs);
%
% cPar.Method='HC';
% cPar.n=6;
%
% cPar.RigidConstraints=indFix(indKeep);
% % cPar.RigidConstraints=cPar.RigidConstraints(cPar.RigidConstraints>0);
%
% [Vb_clean]=tesSmooth(Fb_clean,Vb_clean,[],cPar);
% ind=Fb(:);
% ind=unique(ind(:));
% Vs(ind,:)=Vb_clean;

% cFigure; hold on;
% gpatch(Fb,Vs,'bw','k',1);
% % patchNormPlot(Fs,Vs);
% % plotV(Vs(indKeep,:),'k.','MarkerSize',25)
% axisGeom(gca,fontSize);
% camlight headlight; lighting flat;
% drawnow;

Visualizing input mesh and lattic structures

cFigure;
hs=subplot(1,2,1);
title('The input mesh','fontSize',fontSize)
hold on;
gpatch(Fr,Vr,0.5*ones(1,3),'k',0.5);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

subplot(1,2,2);
title('Lattice side 1','fontSize',fontSize)
hold on;
gpatch(Fb,V,'bw');
% patchNormPlot(Fs,Vs);
axisGeom(gca,fontSize);
camlight headlight; lighting flat;

drawnow;

Visualize BC's

cFigure; hold on;
title('Boundary conditions','FontSize',fontSize);
gpatch(Fb,V,'kw','none',0.4);
hl=gobjects(1,6);
plotColors=gjet(6);
for q=1:1:numel(bcPrescribeListCell)
    hl(q)=plotV(V(bcPrescribeListCell{q},:),'k.','MarkerSize',markerSize);
    hl(q).Color=plotColors(q,:);
end

legend(hl,{'BC 1','BC 2','BC 3','BC 4','BC 5','BC 6'});
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

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='3.0';

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

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
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;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
% febio_spec.Material.material{1}.c2=c1;
% febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;


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

% -> Elements
partName1='Part1_lattice';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> NodeSets
for q=1:1:numel(bcPrescribeListCell)
    febio_spec.Mesh.NodeSet{q}.ATTR.name=['bcPrescribeList_',num2str(q)];
    febio_spec.Mesh.NodeSet{q}.node.ATTR.id=bcPrescribeListCell{q}';
end

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section

% -> Prescribe boundary conditions
directionStringSet={'x','x','y','y','z','z'};
displacementMagnitudeDir=[1 -1 1 -1 1 -1];
for q=1:1:numel(bcPrescribeListCell)
    nodeSetName=febio_spec.Mesh.NodeSet{q}.ATTR.name;
    febio_spec.Boundary.bc{q}.ATTR.type='prescribe';
    febio_spec.Boundary.bc{q}.ATTR.node_set=nodeSetName;
    febio_spec.Boundary.bc{q}.dof=directionStringSet{q};
    febio_spec.Boundary.bc{q}.scale.ATTR.lc=1;
    febio_spec.Boundary.bc{q}.scale.VAL=displacementMagnitudeDir(q).*displacementMagnitude;
    febio_spec.Boundary.bc{q}.relative=0;
end

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%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{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

Quick viewing of the FEBio input file structure

The febView function can be used to view the xml structure in a MATLAB figure window.

febView(febio_spec); %Viewing the febio file

Exporting the FEBio 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
%system(['gedit ',febioFebFileName,' &']);

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.runMode=runMode;
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-------->    RUNNING/MONITORING FEBIO JOB    <-------- 09-Sep-2021 17:05:28
FEBio path: /home/kevin/FEBioStudio/bin/febio3
# Attempt removal of existing log files                09-Sep-2021 17:05:28
 * Removal succesful                                   09-Sep-2021 17:05:28
# Attempt removal of existing .xplt files              09-Sep-2021 17:05:28
 * Removal succesful                                   09-Sep-2021 17:05:28
# Starting FEBio...                                    09-Sep-2021 17:05:28
  Max. total analysis time is: Inf s
 * Waiting for log file creation                       09-Sep-2021 17:05:28
   Max. wait time: 10 s
 * Log file found.                                     09-Sep-2021 17:05:29
# Parsing log file...                                  09-Sep-2021 17:05:29
    number of iterations   : 2                         09-Sep-2021 17:05:30
    number of reformations : 2                         09-Sep-2021 17:05:30
------- converged at time : 0.02                       09-Sep-2021 17:05:30
    number of iterations   : 2                         09-Sep-2021 17:05:31
    number of reformations : 2                         09-Sep-2021 17:05:31
------- converged at time : 0.04                       09-Sep-2021 17:05:31
    number of iterations   : 2                         09-Sep-2021 17:05:33
    number of reformations : 2                         09-Sep-2021 17:05:33
------- converged at time : 0.06                       09-Sep-2021 17:05:33
    number of iterations   : 31                        09-Sep-2021 17:05:48
    number of reformations : 31                        09-Sep-2021 17:05:48
------- converged at time : 0.08                       09-Sep-2021 17:05:48
    number of iterations   : 12                        09-Sep-2021 17:05:57
    number of reformations : 12                        09-Sep-2021 17:05:57
------- converged at time : 0.0949841                  09-Sep-2021 17:05:57
    number of iterations   : 12                        09-Sep-2021 17:06:04
    number of reformations : 12                        09-Sep-2021 17:06:04
------- converged at time : 0.110971                   09-Sep-2021 17:06:04
    number of iterations   : 16                        09-Sep-2021 17:06:12
    number of reformations : 16                        09-Sep-2021 17:06:12
------- converged at time : 0.127761                   09-Sep-2021 17:06:12
    number of iterations   : 5                         09-Sep-2021 17:06:15
    number of reformations : 5                         09-Sep-2021 17:06:15
------- converged at time : 0.145193                   09-Sep-2021 17:06:15
    number of iterations   : 4                         09-Sep-2021 17:06:17
    number of reformations : 4                         09-Sep-2021 17:06:17
------- converged at time : 0.163139                   09-Sep-2021 17:06:17
    number of iterations   : 4                         09-Sep-2021 17:06:19
    number of reformations : 4                         09-Sep-2021 17:06:19
------- converged at time : 0.181495                   09-Sep-2021 17:06:19
    number of iterations   : 4                         09-Sep-2021 17:06:21
    number of reformations : 4                         09-Sep-2021 17:06:21
------- converged at time : 0.20018                    09-Sep-2021 17:06:21
    number of iterations   : 4                         09-Sep-2021 17:06:24
    number of reformations : 4                         09-Sep-2021 17:06:24
------- converged at time : 0.219128                   09-Sep-2021 17:06:24
    number of iterations   : 7                         09-Sep-2021 17:06:28
    number of reformations : 7                         09-Sep-2021 17:06:28
------- converged at time : 0.238287                   09-Sep-2021 17:06:28
    number of iterations   : 5                         09-Sep-2021 17:06:31
    number of reformations : 5                         09-Sep-2021 17:06:31
------- converged at time : 0.257613                   09-Sep-2021 17:06:31
    number of iterations   : 4                         09-Sep-2021 17:06:33
    number of reformations : 4                         09-Sep-2021 17:06:33
------- converged at time : 0.277075                   09-Sep-2021 17:06:33
    number of iterations   : 6                         09-Sep-2021 17:06:36
    number of reformations : 6                         09-Sep-2021 17:06:36
------- converged at time : 0.296644                   09-Sep-2021 17:06:36
    number of iterations   : 9                         09-Sep-2021 17:06:41
    number of reformations : 9                         09-Sep-2021 17:06:41
------- converged at time : 0.316299                   09-Sep-2021 17:06:41
    number of iterations   : 5                         09-Sep-2021 17:06:44
    number of reformations : 5                         09-Sep-2021 17:06:44
------- converged at time : 0.336024                   09-Sep-2021 17:06:44
    number of iterations   : 4                         09-Sep-2021 17:06:46
    number of reformations : 4                         09-Sep-2021 17:06:46
------- converged at time : 0.355803                   09-Sep-2021 17:06:46
    number of iterations   : 4                         09-Sep-2021 17:06:50
    number of reformations : 4                         09-Sep-2021 17:06:50
------- converged at time : 0.375626                   09-Sep-2021 17:06:50
    number of iterations   : 4                         09-Sep-2021 17:06:52
    number of reformations : 4                         09-Sep-2021 17:06:52
------- converged at time : 0.395485                   09-Sep-2021 17:06:52
    number of iterations   : 4                         09-Sep-2021 17:06:55
    number of reformations : 4                         09-Sep-2021 17:06:55
------- converged at time : 0.415372                   09-Sep-2021 17:06:55
    number of iterations   : 4                         09-Sep-2021 17:06:57
    number of reformations : 4                         09-Sep-2021 17:06:57
------- converged at time : 0.435282                   09-Sep-2021 17:06:57
    number of iterations   : 4                         09-Sep-2021 17:07:00
    number of reformations : 4                         09-Sep-2021 17:07:00
------- converged at time : 0.45521                    09-Sep-2021 17:07:00
    number of iterations   : 4                         09-Sep-2021 17:07:02
    number of reformations : 4                         09-Sep-2021 17:07:02
------- converged at time : 0.475152                   09-Sep-2021 17:07:02
    number of iterations   : 4                         09-Sep-2021 17:07:05
    number of reformations : 4                         09-Sep-2021 17:07:05
------- converged at time : 0.495106                   09-Sep-2021 17:07:05
    number of iterations   : 4                         09-Sep-2021 17:07:07
    number of reformations : 4                         09-Sep-2021 17:07:07
------- converged at time : 0.515069                   09-Sep-2021 17:07:07
    number of iterations   : 3                         09-Sep-2021 17:07:09
    number of reformations : 3                         09-Sep-2021 17:07:09
------- converged at time : 0.535039                   09-Sep-2021 17:07:09
    number of iterations   : 3                         09-Sep-2021 17:07:11
    number of reformations : 3                         09-Sep-2021 17:07:11
------- converged at time : 0.555015                   09-Sep-2021 17:07:11
    number of iterations   : 3                         09-Sep-2021 17:07:13
    number of reformations : 3                         09-Sep-2021 17:07:13
------- converged at time : 0.574996                   09-Sep-2021 17:07:13
    number of iterations   : 3                         09-Sep-2021 17:07:14
    number of reformations : 3                         09-Sep-2021 17:07:14
------- converged at time : 0.594981                   09-Sep-2021 17:07:14
    number of iterations   : 3                         09-Sep-2021 17:07:16
    number of reformations : 3                         09-Sep-2021 17:07:16
------- converged at time : 0.614969                   09-Sep-2021 17:07:16
    number of iterations   : 3                         09-Sep-2021 17:07:18
    number of reformations : 3                         09-Sep-2021 17:07:18
------- converged at time : 0.634959                   09-Sep-2021 17:07:18
    number of iterations   : 3                         09-Sep-2021 17:07:20
    number of reformations : 3                         09-Sep-2021 17:07:20
------- converged at time : 0.654952                   09-Sep-2021 17:07:20
    number of iterations   : 3                         09-Sep-2021 17:07:22
    number of reformations : 3                         09-Sep-2021 17:07:22
------- converged at time : 0.674945                   09-Sep-2021 17:07:22
    number of iterations   : 3                         09-Sep-2021 17:07:24
    number of reformations : 3                         09-Sep-2021 17:07:24
------- converged at time : 0.69494                    09-Sep-2021 17:07:24
    number of iterations   : 3                         09-Sep-2021 17:07:26
    number of reformations : 3                         09-Sep-2021 17:07:26
------- converged at time : 0.714936                   09-Sep-2021 17:07:26
    number of iterations   : 3                         09-Sep-2021 17:07:28
    number of reformations : 3                         09-Sep-2021 17:07:28
------- converged at time : 0.734933                   09-Sep-2021 17:07:28
    number of iterations   : 3                         09-Sep-2021 17:07:30
    number of reformations : 3                         09-Sep-2021 17:07:30
------- converged at time : 0.754931                   09-Sep-2021 17:07:30
    number of iterations   : 3                         09-Sep-2021 17:07:32
    number of reformations : 3                         09-Sep-2021 17:07:32
------- converged at time : 0.774929                   09-Sep-2021 17:07:32
    number of iterations   : 3                         09-Sep-2021 17:07:34
    number of reformations : 3                         09-Sep-2021 17:07:34
------- converged at time : 0.794927                   09-Sep-2021 17:07:34
    number of iterations   : 3                         09-Sep-2021 17:07:36
    number of reformations : 3                         09-Sep-2021 17:07:36
------- converged at time : 0.814926                   09-Sep-2021 17:07:36
    number of iterations   : 3                         09-Sep-2021 17:07:38
    number of reformations : 3                         09-Sep-2021 17:07:38
------- converged at time : 0.834925                   09-Sep-2021 17:07:38
    number of iterations   : 3                         09-Sep-2021 17:07:40
    number of reformations : 3                         09-Sep-2021 17:07:40
------- converged at time : 0.854924                   09-Sep-2021 17:07:40
    number of iterations   : 3                         09-Sep-2021 17:07:42
    number of reformations : 3                         09-Sep-2021 17:07:42
------- converged at time : 0.874923                   09-Sep-2021 17:07:42
    number of iterations   : 3                         09-Sep-2021 17:07:44
    number of reformations : 3                         09-Sep-2021 17:07:44
------- converged at time : 0.894923                   09-Sep-2021 17:07:44
    number of iterations   : 3                         09-Sep-2021 17:07:46
    number of reformations : 3                         09-Sep-2021 17:07:46
------- converged at time : 0.914922                   09-Sep-2021 17:07:46
    number of iterations   : 3                         09-Sep-2021 17:07:48
    number of reformations : 3                         09-Sep-2021 17:07:48
------- converged at time : 0.934922                   09-Sep-2021 17:07:48
    number of iterations   : 3                         09-Sep-2021 17:07:50
    number of reformations : 3                         09-Sep-2021 17:07:50
------- converged at time : 0.954922                   09-Sep-2021 17:07:50
    number of iterations   : 3                         09-Sep-2021 17:07:51
    number of reformations : 3                         09-Sep-2021 17:07:51
------- converged at time : 0.974921                   09-Sep-2021 17:07:51
    number of iterations   : 3                         09-Sep-2021 17:07:54
    number of reformations : 3                         09-Sep-2021 17:07:54
------- converged at time : 0.994921                   09-Sep-2021 17:07:54
    number of iterations   : 2                         09-Sep-2021 17:07:56
    number of reformations : 2                         09-Sep-2021 17:07:56
------- converged at time : 1                          09-Sep-2021 17:07:56
 Elapsed time : 0:02:26                                09-Sep-2021 17:07:56
 N O R M A L   T E R M I N A T I O N
# Done                                                 09-Sep-2021 17:07:56
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Import FEBio results

if runFlag==1 %i.e. a succesful run
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
    % Importing nodal forces
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_force),1,1);
    N_force_mat=dataStruct.data;

    indicesSide=bcPrescribeListCell{1};
    areaSide=sampleSize.^2;

    stressVal=mean(squeeze(N_force_mat(indicesSide,1,:))./areaSide,1);
    J_Val=1-((1-J_final).*timeVec(:));
    cFigure; hold on;
    xlabel('$J$','Interpreter','Latex'); ylabel('$\sigma^*$','Interpreter','Latex');
    plot(J_Val(:),stressVal(:),'r.-','LineWidth',3,'MarkerSize',15);
    axis square; axis tight; grid on; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

Plotting the simulated results using anim8 to visualize and animate deformations

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    %     gpatch(Fb,Vs,'kw','none',0.25); %A static graphics object
    hp.FaceColor='interp';

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    %     view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;
end

GIBBON www.gibboncode.org

Kevin Mattheus 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-2021 Kevin Mattheus Moerman and the GIBBON contributors

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