DEMO_febio_0050_foot_insole_01

Below is a demonstration for:

Contents

Keywords

clear; close all; clc;

Plot settings

fontSize=15;
faceAlpha1=1;
faceAlpha2=0.3;
markerSize1=15;
markerSize2=10;
lineWidth=2;

Control parameters

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

% 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

% Surface model file names
fileNameBones=fullfile(loadPathSurfaces,'Foot_bulk.stl');
fileNameSkin=fullfile(loadPathSurfaces,'Skin_coarse.stl');

%Geometric parameters
soleOffsetOutward=3; %How much the insole protrudes outward from the foot
numSmoothIterations_sole_Z=15; %Number of smoothing iterations for the sole Z value
soleMinThickness=6; %Sole thickness
volumeFactor=1; %Volume factor used in tetgen, larger means larger internal elements
footSize=251; %Foot size in mm size 6=241 mm, size 7=251 mm

% Material parameters
%-> Tissue (Ogden, uncoupled)
% Data source: Petre et al. 2013, Optimization of Nonlinear Hyperelastic
% Coefficients for Foot Tissues Using a Magnetic Resonance Imaging
% Deformation Experiment. https://dx.doi.org/10.1115%2F1.4023695
c1_1=0.02652; %Shear-modulus-like parameter in MPa
m1_1=17.71; %Material parameter setting degree of non-linearity
k_1=c1_1*100; %Bulk modulus

%-> Sole
%Sole material parameters (Ogden unconstrained, with m=2 i.e. Neo-Hookean)
c1_min=0.05; %Shear-modulus-like parameter in MPa
c1_max=0.1; %Shear-modulus-like parameter in MPa
numMaterialLevels=50; %Max number of materials to use
initialMaterialLevel=numMaterialLevels;
c1_range=linspace(c1_min,c1_max,numMaterialLevels); %Shear-modulus-like parameter
bulkModulusFactor=10;

testCase=1; % Change to test different sole cases, e.g. stiff, soft, spatially varying.
maxLevelColorbar_SED=0.0005;

% 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='internal';

initialSoleSpacing=0.1;

%Boundary condition parameters
displacementMagnitude=-2.63; %Displacement applied to the bones
bodyWeight=65/2;% Body weight in Kg
forceBody=bodyWeight.*-9.81; %Force due to body weight
forceDifferenceToleranceFraction=0.05;
forceDifferenceTolerance=abs(forceBody.*forceDifferenceToleranceFraction);
optimizeForceOption=1; %set to 0 to turn off and use input displacement
maxProposedDisplacementDifference=1;

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

Import surface models

%Import STL file for the bones
[stlStruct] = import_STL(fileNameBones);
F1=stlStruct.solidFaces{1}; %Faces
V1=stlStruct.solidVertices{1}; %Vertices
[F1,V1]=mergeVertices(F1,V1); %Merge nodes

%Import STL file for the skin
[stlStruct] = import_STL(fileNameSkin);
F2=stlStruct.solidFaces{1}; %Faces
V2=stlStruct.solidVertices{1}; %Vertices
[F2,V2]=mergeVertices(F2,V2); %Merge nodes
V2=V2*100; %Scale

Subdivide skin mesh and smooth

% %Split elements once
% [F2,V2]=subtri(F2,V2,1);
%
% %Smooth
% clear cParSmooth;
% cParSmooth.n=5;
% cParSmooth.Method='HC';
% [V2]=patchSmooth(F2,V2,[],cParSmooth);

Visualize imported surfaces

cFigure; hold on;
gpatch(F1,V1,'bw','none',1);
gpatch(F2,V2,'rw','none',0.5);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Compute mesh derived parameters

Compute point spacings for the surfaces meshes. These are useful to relate other mesh sizes to.

pointSpacing1=mean(patchEdgeLengths(F1,V1));
pointSpacing2=mean(patchEdgeLengths(F2,V2));
pointSpacing=mean([pointSpacing1 pointSpacing2]);
pointSpacingSole=pointSpacing;

Reorient surfaces

Reorient so that the leg points up in a forward view with the toes pointing towards the viewer.

R=euler2DCM([-0.5*pi 0 -0.5*pi]);
V1=V1*R;
V2=V2*R;

Visualize reoriented surfaces

cFigure; hold on;
gpatch(F1,V1,'bw','none',1);
gpatch(F2,V2,'rw','none',0.5);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Cut skin surface

The skin surface is cut such that it stops where the bones of the ankle end in the z direction.

%Create a logic for cutting away faces
max_Z1=max(V1(:,3))+2*pointSpacing; %Max Z-level used for cutting
logicVertices=V2(:,3)<max_Z1; %Logic for the points below this level
logicFaces=all(logicVertices(F2),2); %Logic for the faces
logicFaces=triSurfLogicSharpFix(F2,logicFaces,3); %Altered logic so it is smoother

Visualize

cFigure; hold on;
gpatch(F1,V1,'bw','none',1);
gpatch(F2,V2,logicFaces,'none',0.5);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Cut away faces using logic

F2=F2(logicFaces,:); %The faces to keep
[F2,V2]=patchCleanUnused(F2,V2); %Remove unused points

%Attempt to self triangulate potentially jagged edge
Eb=patchBoundary(F2,V2); %Get boundary edges
indBoundary=edgeListToCurve(Eb); %Convert boundary edges to a curve list
indBoundary=indBoundary(1:end-1); %Trim off last point since it is equal to first on a closed loop
angleThreshold=pi*(120/180); %threshold for self triangulation
[F2,V2,indBoundaryTop]=triSurfSelfTriangulateBoundary(F2,V2,indBoundary,angleThreshold,1);

%Force boundary to have the max Z level chosen
V2(indBoundaryTop,3)=max_Z1;

Visualize

cFigure; hold on;
gpatch(F1,V1,'bw','none',1);
gpatch(F2,V2,'rw','k',0.5);
plotV(V2(indBoundaryTop,:),'r-','LineWidth',lineWidth);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Scale foot to desired size

currentFootSize=abs(max(V2(:,1))-min(V2(:,1)));
V2=V2./currentFootSize; %Scale to unit foot
V2=V2.*footSize; %Scale to desired foot size
V1=V1./currentFootSize; %Scale to unit foot
V1=V1.*footSize; %Scale to desired foot size

Close over top of skin

The top boundary curve of the cut surface is filled with triangles. This is a 2D method. The z-coordinate is added after.

[F2t,V2t]=regionTriMesh2D({V2(indBoundaryTop,[1 2])},pointSpacing2,0,0);
V2t(:,3)=mean(V2(indBoundaryTop,3)); %Add/set z-level

Visualize

cFigure; hold on;
gpatch(F1,V1,'bw','none',1);
gpatch(F2,V2,'rw','k',0.5);
gpatch(F2t,V2t,'gw','k',1);
plotV(V2(indBoundaryTop,:),'r-','LineWidth',lineWidth);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Joining surface features

Add all surface sets together in joint list of faces, vertices

[FT,VT,CT]=joinElementSets({F1,F2,F2t},{V1,V2,V2t});

Merge shared nodes

The join operation only adds the sets together. Nodes with the same coordinates are not seen as the same yet and need to be merged.

[FT,VT]=mergeVertices(FT,VT); %Merge nodes

Visualize

cFigure; hold on;
gpatch(FT,VT,CT,'k',0.5);
patchNormPlot(FT,VT);
camlight('headlight');
axisGeom(gca,fontSize);
colormap gjet; icolorbar;
drawnow;

Mesh foot 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({FT(CT==2 | CT==3,:),FT(CT==1,:)},{VT,VT});
[V_hole]=getInnerPoint(FT(CT==1,:),VT);

Visualize interior points

cFigure; hold on;
gpatch(FT,VT,'kw','none',0.2);
hp1=plotV(V_region,'r.','markerSize',markerSize1);
hp2=plotV(V_hole,'b.','markerSize',markerSize1);
legend([hp1 hp2],{'Region point','Hole point'});
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

Mesh using tetgen

inputStruct.stringOpt='-pq1.2AaY'; %TetGen option string
inputStruct.Faces=FT; %The faces
inputStruct.Nodes=VT; %The vertices
inputStruct.holePoints=V_hole; %The hole interior points
inputStruct.faceBoundaryMarker=CT; %Face boundary markers
inputStruct.regionPoints=V_region; %The region interior points
inputStruct.regionA=tetVolMeanEst(F2,V2)*volumeFactor; %Volume for regular tets

Mesh model using tetrahedral elements using tetGen

[meshOutput]=runTetGen(inputStruct); %Run tetGen
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- TETGEN Tetrahedral meshing --- 10-Aug-2020 10:37:29
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing SMESH file --- 10-Aug-2020 10:37:29
----> Adding node field
----> Adding facet field
----> Adding holes specification
----> Adding region specification
--- Done --- 10-Aug-2020 10:37:29
--- Running TetGen to mesh input boundary--- 10-Aug-2020 10:37:29
Opening /mnt/data/MATLAB/GIBBON/data/temp/temp.smesh.
Delaunizing vertices...
Delaunay seconds:  0.038629
Creating surface mesh ...
Surface mesh seconds:  0.012508
Recovering boundaries...
Boundary recovery seconds:  0.031729
Removing exterior tetrahedra ...
Spreading region attributes.
Exterior tets removal seconds:  0.015522
Recovering Delaunayness...
Delaunay recovery seconds:  0.011733
Refining mesh...
Refinement seconds:  0.154201
Optimizing mesh...
Optimization seconds:  0.008934

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.142473
Total running seconds:  0.416277

Statistics:

  Input points: 6538
  Input facets: 13068
  Input segments: 19602
  Input holes: 1
  Input regions: 1

  Mesh points: 10942
  Mesh tetrahedra: 48791
  Mesh faces: 104116
  Mesh faces on exterior boundary: 13068
  Mesh faces on input facets: 13068
  Mesh edges on input segments: 19602
  Steiner points inside domain: 4404

--- Done --- 10-Aug-2020 10:37:30
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Importing TetGen files --- 10-Aug-2020 10:37:30
--- Done --- 10-Aug-2020 10:37:30

Access model element and patch data

Fb_foot=meshOutput.facesBoundary; %Boundary faces of the foot
Cb_foot=meshOutput.boundaryMarker; %Boundary marker/color data for the foot
V_foot=meshOutput.nodes; %The vertices/nodes
E_foot=meshOutput.elements; %The tet4 elements

Visualizing mesh using meshView, see also anim8

meshView(meshOutput);

Build insole

The insole is created by taking the 2D convex hull of the foot (ignoring the z-direction). Next this convex hull is resampled and filled with triangular elements. The z-coordinates are then based on the nearest foot nodes. The z-coordinate data is next smoothed to create a smooth surface. This surface is next thickened to form hexahedral elements.

Create sole boundary curve in 2D

% Offset surface outward to thicken so sole is enlarged outward
[~,~,Nv]=patchNormal(FT,VT);
VT_sole=VT+soleOffsetOutward.*Nv;

P=VT_sole(:,[1 2]); %Point set flattened to 2D
DT=delaunayTriangulation(P); %Delaunay triangulation of 2D set
VD=DT.Points; %Delaunay point set
VD(:,3)=min(VT(:,3)); %Set z-coord to minimum for now
indChull=DT.convexHull; %Ordered point list for conhex hull
indChull=indChull(1:end-1); %Trim away last (=start) point to avoid double
V_chull=VD(indChull,:); %Vertices for convex hull
D=max(pathLength(V_chull)); %Get length of sole curve for resampling
numResample=ceil(D./pointSpacingSole);
V_sole_curve=evenlySampleCurve(V_chull,numResample,'pchip',1);

Visualize

cFigure; hold on;
gpatch(FT,VT,'kw','none',0.25);
plotV(P,'k.','MarkerSize',markerSize2);
plotV(V_chull,'r.','MarkerSize',markerSize1*2);
plotV(V_sole_curve,'b.-','MarkerSize',markerSize1,'LineWidth',lineWidth);
camlight('headlight');
axisGeom(gca,fontSize);
colormap gjet; icolorbar;
drawnow;

Build sole top surface

[F_sole_top,V_sole_top]=regionTriMesh2D({V_sole_curve(:,[1 2])},pointSpacingSole,0,0);
V_sole_top(:,3)=mean(V_sole_curve(:,3));

Eb=patchBoundary(F_sole_top,V_sole_top);
indBoundary=unique(Eb(:));

%Get z-coordinate
[~,indMin]=minDist(V_sole_top,VT);
V_sole_top(:,3)=VT(indMin,3);

%Free smoothing of boundary
clear cParSmooth;
cParSmooth.n=numSmoothIterations_sole_Z;
cParSmooth.Method='LAP';
[p]=patchSmooth(Eb,V_sole_top,[],cParSmooth);
V_sole_top(indBoundary,3)=p(indBoundary,3);

% Constrained general smoothing
cParSmooth.n=numSmoothIterations_sole_Z;
cParSmooth.Method='LAP';
cParSmooth.RigidConstraints=indBoundary;
[V_sole_top(:,3)]=patchSmooth(F_sole_top,V_sole_top(:,3),[],cParSmooth);

%Shift a small amount to there is no initial contact
[~,indMin]=minDist(V_sole_top,VT);
dz=VT(indMin,3)-V_sole_top(:,3);
V_sole_top(:,3)=V_sole_top(:,3)-abs(min(dz));

V_sole_top(:,3)=V_sole_top(:,3)-initialSoleSpacing;

Visualize

cFigure; hold on;
gpatch(FT,VT,'kw','none',0.25);
gpatch(F_sole_top,V_sole_top,'gw','k',1);
patchNormPlot(F_sole_top,V_sole_top);
camlight('headlight');
axisGeom(gca,fontSize);
colormap gjet; icolorbar;
drawnow;

Create bottom surface of sole

dirSet=[0 0 -1];
numElementsSoleThickness=ceil(soleMinThickness./pointSpacingSole);
[E_sole,V_sole,F_sole_top,F_sole_bottom]=patchThick(fliplr(F_sole_top),V_sole_top,dirSet,soleMinThickness,numElementsSoleThickness);
F_sole_top=fliplr(F_sole_top);

% Use element2patch to get patch data
F_E_sole=element2patch(E_sole,[],'penta6');
% indBoundaryFaces=tesBoundary(F_E_sole,V_sole);
% Fb_sole=F_E_sole(indBoundaryFaces,:);
cFigure; hold on;
gpatch(FT,VT,'kw','none',0.25);
gpatch(F_sole_top,V_sole,'gw','k',1);
gpatch(F_sole_bottom,V_sole,'bw','k',1);
patchNormPlot(F_sole_top,V_sole);
camlight('headlight');
axisGeom(gca,fontSize);
colormap gjet; icolorbar;
drawnow;

Visualize

cFigure; hold on;
title('Hexahedral mesh');
gpatch(FT,VT,'kw','none',0.25);
gpatch(F_E_sole,V_sole,'bw','k',1);
% patchNormPlot(Fb_sole,V_sole);
axisGeom;
camlight headlight;
drawnow;

Joining node sets

V=[V_foot;V_sole;]; %Combined node sets
E_sole=E_sole+size(V_foot,1); %Fixed element indices
F_sole_top=F_sole_top+size(V_foot,1); %Fixed element indices
F_sole_bottom=F_sole_bottom+size(V_foot,1); %Fixed indices
F_E_sole=element2patch(E_sole,[],'penta6');

Visualize

cFigure; hold on;
title('Hexahedral mesh');
gpatch(Fb_foot,V,Cb_foot,'k',1);
gpatch(F_sole_top,V,'kw','k',1);
gpatch(F_sole_bottom,V,'kw','k',1);
% patchNormPlot(FEs,V);
colormap gjet; icolorbar;
axisGeom;
camlight headlight;
drawnow;

Define contact surfaces

% The rigid master surface of the sphere
F_contact_master=F_sole_top;

% The deformable slave surface of the slab
F_contact_slave=fliplr(Fb_foot(Cb_foot==2,:));

Visualize contact surfaces

cFigure; hold on;
title('Contact sets and normal directions','FontSize',fontSize);

gpatch(Fb_foot,V,'kw','none',faceAlpha2);
% gpatch(Fb_sole,V,'kw','none',faceAlpha2);

hl(1)=gpatch(F_contact_master,V,'gw','k',1);
patchNormPlot(F_contact_master,V);
hl(2)=gpatch(F_contact_slave,V,'bw','k',1);
patchNormPlot(F_contact_slave,V);

legend(hl,{'Master','Slave'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

Define boundary conditions

%Supported nodes
bcSupportList=unique(F_sole_bottom);

%Prescribed displacement nodes
bcPrescribeList=unique(Fb_foot(Cb_foot==1,:));

Visualize BC's

hf=cFigure; hold on;
title('Boundary conditions model','FontSize',fontSize);
gpatch(Fb_foot,V,'kw','none',faceAlpha2);

hl2(1)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize2);
hl2(2)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize2);
legend(hl2,{'BC prescribe','BC support'});
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

Visualize rigid body bone surface elements

F_bone=Fb_foot(Cb_foot==1,:);

center_of_mass_bone=mean(V(bcPrescribeList,:),1);

cFigure; hold on;
% title('Boundary conditions model','FontSize',fontSize);
gpatch(Fb_foot,V,'kw','none',faceAlpha2);
gpatch(F_bone,V,'w','k',1);
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

Set-up materials

switch testCase
    case 0 %Spatially varying based on bone distance
        VE=patchCentre(E_sole,V); %Element centres
        V_vulnerable=VT(bcPrescribeList,:);
        W=minDist(VE,V_vulnerable); %Distance to bone points
        W=W-min(W(:));
        W=W./max(W(:));
        W=W.*(c1_max-c1_min);
        W=W+c1_min;
        c1_desired=W;
    case 1 %
        c1_desired=c1_max*ones(size(E_sole,1),1);
    case 2
        c1_desired=(c1_max+c1_min)/2*ones(size(E_sole,1),1);
    case 3
        c1_desired=c1_min*ones(size(E_sole,1),1);
end

%Snap to nearest available materials
[~,elementMaterialID]=minDist(c1_desired(:),c1_range(:));

Visualize material levels

%Use element2patch to get patch and color data
[F_E_sole,C_F_E_sole]=element2patch(E_sole,elementMaterialID,'penta6');

hf=cFigure; hold on;
title('Material levels','FontSize',fontSize);
gpatch(Fb_foot(Cb_foot==1,:),V,'kw','none',1);
gpatch(Fb_foot(Cb_foot~=1,:),V,'w','none',0.5);


for q=1:1:numel(F_E_sole)
    gpatch(F_E_sole{q},V,C_F_E_sole{q},'k',1);
end
axisGeom(gca,fontSize);
camlight headlight;
colormap parula;
icolorbar;
%caxis([0 maxLevelColorbar_SED]);
drawnow;

Convert list of materials to element id

%Sorting elements according to material label
[elementMaterialID,indSortElements]=sort(elementMaterialID);
E_sole=E_sole(indSortElements,:);

%Removing unused materials
[indUsed,ind1,ind2]=unique(elementMaterialID(:));
elementMaterialID=ind2;
c1=c1_range(indUsed);

numMaterials=numel(indUsed);

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}.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=center_of_mass_bone;

%Sole materials
for q=1:1:numMaterials
    febio_spec.Material.material{q+2}.ATTR.type='Ogden unconstrained';
    febio_spec.Material.material{q+2}.ATTR.id=q+2;
    febio_spec.Material.material{q+2}.c1=c1(q);
    febio_spec.Material.material{q+2}.m1=2; %m1=2 to make it Neo-Hookean
    febio_spec.Material.material{q+2}.cp=c1(q)*bulkModulusFactor;
end

%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='Foot'; %Name of the element set
febio_spec.Geometry.Elements{1}.elem.ATTR.id=(1:1:size(E_foot,1))'; %Element id's
febio_spec.Geometry.Elements{1}.elem.VAL=E_foot;

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='Bone'; %Name of the element set
febio_spec.Geometry.Elements{2}.elem.ATTR.id=(size(E_foot,1)+1:1:size(E_foot,1)+size(F_bone,1))'; %Element id's
febio_spec.Geometry.Elements{2}.elem.VAL=F_bone;

n=size(E_foot,1)+size(F_bone,1)+1;
for q=1:1:numMaterials
    logicMaterialNow=(elementMaterialID==q);
    febio_spec.Geometry.Elements{q+2}.ATTR.type='penta6'; %Element type of this set
    febio_spec.Geometry.Elements{q+2}.ATTR.mat=q+2; %material index for this set
    febio_spec.Geometry.Elements{q+2}.ATTR.name=['Sole_mat',num2str(q)]; %Name of the element set
    febio_spec.Geometry.Elements{q+2}.elem.ATTR.id=(n:1:(n-1+nnz(logicMaterialNow)))'; %Element id's
    febio_spec.Geometry.Elements{q+2}.elem.VAL=E_sole(logicMaterialNow,:);
    n=n+nnz(logicMaterialNow);
end

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

febio_spec.Geometry.NodeSet{2}.ATTR.name='bcPrescribeList';
febio_spec.Geometry.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

% -> Surfaces
febio_spec.Geometry.Surface{1}.ATTR.name='contact_master';
febio_spec.Geometry.Surface{1}.tri3.ATTR.lid=(1:1:size(F_contact_master,1))';
febio_spec.Geometry.Surface{1}.tri3.VAL=F_contact_master;

febio_spec.Geometry.Surface{2}.ATTR.name='contact_slave';
febio_spec.Geometry.Surface{2}.tri3.ATTR.lid=(1:1:size(F_contact_slave,1))';
febio_spec.Geometry.Surface{2}.tri3.VAL=F_contact_slave;

% -> Surface pairs
febio_spec.Geometry.SurfacePair{1}.ATTR.name='Contact1';
febio_spec.Geometry.SurfacePair{1}.master.ATTR.surface=febio_spec.Geometry.Surface{1}.ATTR.name;
febio_spec.Geometry.SurfacePair{1}.slave.ATTR.surface=febio_spec.Geometry.Surface{2}.ATTR.name;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.fix{1}.ATTR.bc='x';
febio_spec.Boundary.fix{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Boundary.fix{2}.ATTR.bc='y';
febio_spec.Boundary.fix{2}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Boundary.fix{3}.ATTR.bc='z';
febio_spec.Boundary.fix{3}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;

febio_spec.Boundary.fix{4}.ATTR.bc='x';
febio_spec.Boundary.fix{4}.ATTR.node_set=febio_spec.Geometry.NodeSet{2}.ATTR.name;
febio_spec.Boundary.fix{5}.ATTR.bc='y';
febio_spec.Boundary.fix{5}.ATTR.node_set=febio_spec.Geometry.NodeSet{2}.ATTR.name;

% -> 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='Rx';
febio_spec.Boundary.rigid_body{1}.fixed{4}.ATTR.bc='Ry';
febio_spec.Boundary.rigid_body{1}.fixed{5}.ATTR.bc='Rz';

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

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

%Contact section
febio_spec.Contact.contact{1}.ATTR.surface_pair=febio_spec.Geometry.SurfacePair{1}.ATTR.name;
febio_spec.Contact.contact{1}.ATTR.type='sliding-elastic';
febio_spec.Contact.contact{1}.two_pass=1;
febio_spec.Contact.contact{1}.laugon=laugon;
febio_spec.Contact.contact{1}.tolerance=0.2;
febio_spec.Contact.contact{1}.gaptol=0;
febio_spec.Contact.contact{1}.minaug=minaug;
febio_spec.Contact.contact{1}.maxaug=maxaug;
febio_spec.Contact.contact{1}.search_tol=0.01;
febio_spec.Contact.contact{1}.search_radius=0.1;
febio_spec.Contact.contact{1}.symmetric_stiffness=0;
febio_spec.Contact.contact{1}.auto_penalty=1;
febio_spec.Contact.contact{1}.penalty=contactPenalty;
febio_spec.Contact.contact{1}.fric_coeff=fric_coeff;

%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_foot,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

if optimizeForceOption==0

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!!!!!!!!


    if runFlag~=1
        error('FEBio error');
    end
elseif optimizeForceOption==1

    forceDifference=1e6;
    while 1

Set displacement magnitude in input file structure

        febio_spec.Boundary.rigid_body{1}.prescribed.VAL=displacementMagnitude;

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!!!!!!!!

        if runFlag~=1
            error('FEBio error');
        end
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- STARTING FEBIO JOB --- 10-Aug-2020 10:37:51
===========================================================================
         ________    _________   _________     __     _________            
        |        |\ |        |\ |        |\   |  |\  /         \\          
        |    ____|| |    ____|| |    __  ||   |__|| |    ___    ||         
        |   |\___\| |   |\___\| |   |\_| ||    \_\| |   //  \   ||         
        |   ||      |   ||      |   || | ||    __   |  ||    |  ||         
        |   ||__    |   ||__    |   ||_| ||   |  |\ |  ||    |  ||         
        |       |\  |       |\  |         \\  |  || |  ||    |  ||         
        |    ___||  |    ___||  |    ___   || |  || |  ||    |  ||         
        |   |\__\|  |   |\__\|  |   |\__|  || |  || |  ||    |  ||         
        |   ||      |   ||      |   ||  |  || |  || |  ||    |  ||         
        |   ||      |   ||___   |   ||__|  || |  || |   \\__/   ||         
        |   ||      |        |\ |          || |  || |           ||         
        |___||      |________|| |__________|| |__||  \_________//          
                                                                           
      F I N I T E   E L E M E N T S   F O R   B I O M E C H A N I C S      
                                                                           
                 --- v e r s i o n - 2 . 9 . 1 ---                 
                                                                           
                                                                           
  http://febio.org                                                         
  FEBio is a registered trademark.                                         
  copyright (c) 2006-2019 - All rights reserved                            
                                                                           
===========================================================================


Reading file /mnt/data/MATLAB/GIBBON/data/temp/tempModel.feb ...SUCCESS!
 ]0;(0%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 1 : 0.1 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 998983

 1
 Nonlinear solution status: time= 0.1
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.796895e+04    1.803028e-22    0.000000e+00 
	   energy              1.241280e+03    9.870988e-12    1.241280e+01 
	   displacement        4.110022e+02    4.110022e+02    4.110022e-04 
 *************************************************************************
 *                               WARNING                                 *
 *                                                                       *
 * No force acting on the system.                                        *
 *************************************************************************

------- converged at time : 0.1

 ]0;(10%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 2 : 0.2 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 998983

 1
 Nonlinear solution status: time= 0.2
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.796895e+04    2.086304e+00    0.000000e+00 
	   energy              1.241280e+03    5.566173e-01    1.241280e+01 
	   displacement        4.110022e+02    4.110022e+02    4.110022e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1000315

 2
 Nonlinear solution status: time= 0.2
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.796895e+04    2.406821e-02    0.000000e+00 
	   energy              1.241280e+03    1.025380e-02    1.241280e+01 
	   displacement        4.110022e+02    4.544411e-02    4.102763e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1000099

 3
 Nonlinear solution status: time= 0.2
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.796895e+04    1.909944e-04    0.000000e+00 
	   energy              1.241280e+03    2.823314e-05    1.241280e+01 
	   displacement        4.110022e+02    1.509431e-03    4.101657e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 999883

 4
 Nonlinear solution status: time= 0.2
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.796895e+04    1.935657e-08    0.000000e+00 
	   energy              1.241280e+03    4.461726e-07    1.241280e+01 
	   displacement        4.110022e+02    8.274500e-05    4.101414e-04 

------- converged at time : 0.2

 ]0;(20%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 3 : 0.3 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 999883

 1
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    2.092443e+01    0.000000e+00 
	   energy              1.240982e+03    2.570209e+00    1.240982e+01 
	   displacement        4.081866e+02    4.081866e+02    4.081866e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1003285

 2
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    2.944138e-01    0.000000e+00 
	   energy              1.240982e+03    7.766243e-02    1.240982e+01 
	   displacement        4.081866e+02    1.536039e-01    4.069055e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1002691

 3
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    9.966429e-03    0.000000e+00 
	   energy              1.240982e+03    3.097391e-03    1.240982e+01 
	   displacement        4.081866e+02    3.928415e-02    4.056046e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1002493

 4
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    7.562043e-04    0.000000e+00 
	   energy              1.240982e+03    2.176143e-04    1.240982e+01 
	   displacement        4.081866e+02    2.380103e-03    4.053643e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1002421

 5
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    5.767179e-04    0.000000e+00 
	   energy              1.240982e+03    1.462199e-04    1.240982e+01 
	   displacement        4.081866e+02    7.619656e-04    4.053330e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1002421

 6
 Nonlinear solution status: time= 0.3
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 6
	step from line search         = 0.563210
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.797208e+04    8.071223e-04    0.000000e+00 
	   energy              1.240982e+03    2.979701e-05    1.240982e+01 
	   displacement        4.081866e+02    1.167117e-04    4.053376e-04 

------- converged at time : 0.3

 ]0;(30%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 4 : 0.4 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1002421

 1
 Nonlinear solution status: time= 0.4
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.799120e+04    2.554448e+01    0.000000e+00 
	   energy              1.240559e+03    3.201492e+00    1.240559e+01 
	   displacement        4.021993e+02    4.021993e+02    4.021993e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1006939

 2
 Nonlinear solution status: time= 0.4
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.799120e+04    1.908973e-01    0.000000e+00 
	   energy              1.240559e+03    8.776823e-02    1.240559e+01 
	   displacement        4.021993e+02    3.680662e-01    3.993855e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1005679

 3
 Nonlinear solution status: time= 0.4
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.799120e+04    3.723294e-03    0.000000e+00 
	   energy              1.240559e+03    1.613900e-03    1.240559e+01 
	   displacement        4.021993e+02    5.562036e-02    3.982030e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1005301

 4
 Nonlinear solution status: time= 0.4
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.799120e+04    3.877538e-04    0.000000e+00 
	   energy              1.240559e+03    9.671478e-05    1.240559e+01 
	   displacement        4.021993e+02    1.694844e-03    3.980232e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1005301

 5
 Nonlinear solution status: time= 0.4
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.799120e+04    1.036677e-03    0.000000e+00 
	   energy              1.240559e+03    5.674766e-05    1.240559e+01 
	   displacement        4.021993e+02    9.348651e-05    3.980068e-04 

------- converged at time : 0.4

 ]0;(40%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 5 : 0.5 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1005301

 1
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    1.292483e+02    0.000000e+00 
	   energy              1.240314e+03    9.139384e+00    1.240314e+01 
	   displacement        3.943717e+02    3.943717e+02    3.943717e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1014481

 2
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    8.794516e-01    0.000000e+00 
	   energy              1.240314e+03    2.279597e-01    1.240314e+01 
	   displacement        3.943717e+02    5.148078e-01    3.920744e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1012807

 3
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    2.912116e-02    0.000000e+00 
	   energy              1.240314e+03    8.648126e-03    1.240314e+01 
	   displacement        3.943717e+02    1.177714e-01    3.904469e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011853

 4
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    2.280841e-03    0.000000e+00 
	   energy              1.240314e+03    5.818329e-04    1.240314e+01 
	   displacement        3.943717e+02    8.606176e-03    3.900724e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011547

 5
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    1.400269e-03    0.000000e+00 
	   energy              1.240314e+03    1.322913e-03    1.240314e+01 
	   displacement        3.943717e+02    2.297775e-02    3.897033e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011529

 6
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 6
	step from line search         = 0.608566
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    1.430983e-03    0.000000e+00 
	   energy              1.240314e+03    2.434760e-04    1.240314e+01 
	   displacement        3.943717e+02    8.807878e-03    3.899116e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011529

 7
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 7
	step from line search         = 0.495013
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    1.604957e-05    0.000000e+00 
	   energy              1.240314e+03    2.326681e-05    1.240314e+01 
	   displacement        3.943717e+02    4.240398e-04    3.899283e-04 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011529

 8
 Nonlinear solution status: time= 0.5
	stiffness updates             = 0
	right hand side evaluations   = 11
	stiffness matrix reformations = 8
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.803796e+04    1.088149e-07    0.000000e+00 
	   energy              1.240314e+03    2.342780e-07    1.240314e+01 
	   displacement        3.943717e+02    4.417656e-06    3.899299e-04 

------- converged at time : 0.5

 ]0;(50%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 6 : 0.6 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1011529

 1
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    2.957778e+02    0.000000e+00 
	   energy              1.239935e+03    1.957633e+01    1.239935e+01 
	   displacement        3.866764e+02    3.866764e+02    3.866764e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1031887

 2
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    2.312718e+00    0.000000e+00 
	   energy              1.239935e+03    5.031091e-01    1.239935e+01 
	   displacement        3.866764e+02    1.290585e+00    3.837173e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1028485

 3
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    1.094441e-01    0.000000e+00 
	   energy              1.239935e+03    2.358332e-02    1.239935e+01 
	   displacement        3.866764e+02    2.725838e-01    3.811951e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1026865

 4
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    8.717353e-03    0.000000e+00 
	   energy              1.239935e+03    2.295966e-03    1.239935e+01 
	   displacement        3.866764e+02    4.041140e-02    3.805220e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1025839

 5
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    1.117799e-03    0.000000e+00 
	   energy              1.239935e+03    1.715853e-04    1.239935e+01 
	   displacement        3.866764e+02    8.360944e-03    3.802882e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1025515

 6
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    3.011325e-04    0.000000e+00 
	   energy              1.239935e+03    4.040207e-05    1.239935e+01 
	   displacement        3.866764e+02    7.951177e-04    3.802659e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1025461

 7
 Nonlinear solution status: time= 0.6
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 7
	step from line search         = 0.191079
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.811788e+04    1.971171e-04    0.000000e+00 
	   energy              1.239935e+03    1.566960e-05    1.239935e+01 
	   displacement        3.866764e+02    4.066454e-06    3.802646e-04 

------- converged at time : 0.6

 ]0;(60%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 7 : 0.7 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1025479

 1
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    4.924342e+02    0.000000e+00 
	   energy              1.238636e+03    2.617592e+01    1.238636e+01 
	   displacement        3.735536e+02    3.735536e+02    3.735536e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1047565

 2
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    2.973187e+00    0.000000e+00 
	   energy              1.238636e+03    7.358470e-01    1.238636e+01 
	   displacement        3.735536e+02    1.907894e+00    3.728202e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1043821

 3
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    1.283574e-01    0.000000e+00 
	   energy              1.238636e+03    3.115278e-02    1.238636e+01 
	   displacement        3.735536e+02    3.062244e-01    3.716089e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1042093

 4
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    1.960271e-02    0.000000e+00 
	   energy              1.238636e+03    7.621116e-03    1.238636e+01 
	   displacement        3.735536e+02    1.727029e-01    3.712533e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1040977

 5
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    5.414042e-03    0.000000e+00 
	   energy              1.238636e+03    1.843862e-03    1.238636e+01 
	   displacement        3.735536e+02    7.310127e-02    3.715930e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1040851

 6
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    4.482892e-03    0.000000e+00 
	   energy              1.238636e+03    4.165314e-04    1.238636e+01 
	   displacement        3.735536e+02    1.561368e-02    3.718221e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1040725

 7
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 7
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    1.027832e-03    0.000000e+00 
	   energy              1.238636e+03    6.078621e-04    1.238636e+01 
	   displacement        3.735536e+02    9.023423e-03    3.718190e-04 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1040725

 8
 Nonlinear solution status: time= 0.7
	stiffness updates             = 0
	right hand side evaluations   = 9
	stiffness matrix reformations = 8
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.824840e+04    4.899788e-04    0.000000e+00 
	   energy              1.238636e+03    3.822848e-06    1.238636e+01 
	   displacement        3.735536e+02    2.499099e-04    3.718243e-04 

------- converged at time : 0.7

 ]0;(70%) tempModel.feb - FEBio 2.9.1.0  
===== beginning time step 8 : 0.8 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1040743

 1
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    7.170573e+02    0.000000e+00 
	   energy              1.238064e+03    3.861880e+01    1.238064e+01 
	   displacement        3.719669e+02    3.719669e+02    3.719669e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1064341

 2
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    2.939765e+00    0.000000e+00 
	   energy              1.238064e+03    7.154585e-01    1.238064e+01 
	   displacement        3.719669e+02    2.383462e+00    3.727398e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1060651

 3
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    1.891043e-01    0.000000e+00 
	   energy              1.238064e+03    4.094640e-02    1.238064e+01 
	   displacement        3.719669e+02    4.765248e-01    3.708198e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1058167

 4
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    6.622495e-02    0.000000e+00 
	   energy              1.238064e+03    2.040256e-02    1.238064e+01 
	   displacement        3.719669e+02    4.516836e-01    3.719414e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056745

 5
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    1.975534e-02    0.000000e+00 
	   energy              1.238064e+03    6.911008e-03    1.238064e+01 
	   displacement        3.719669e+02    5.430746e-01    3.748301e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056241

 6
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    8.323247e-03    0.000000e+00 
	   energy              1.238064e+03    7.627506e-06    1.238064e+01 
	   displacement        3.719669e+02    3.616922e-02    3.755985e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056205

 7
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 7
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    3.371058e-03    0.000000e+00 
	   energy              1.238064e+03    2.921891e-04    1.238064e+01 
	   displacement        3.719669e+02    9.156337e-03    3.756741e-04 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056205

 8
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 9
	stiffness matrix reformations = 8
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    3.685015e-03    0.000000e+00 
	   energy              1.238064e+03    1.292385e-04    1.238064e+01 
	   displacement        3.719669e+02    2.324470e-03    3.756752e-04 
Reforming stiffness matrix: reformation #9

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056205

 9
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 9
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    9.965562e-04    0.000000e+00 
	   energy              1.238064e+03    2.565870e-06    1.238064e+01 
	   displacement        3.719669e+02    2.267489e-03    3.756533e-04 
Reforming stiffness matrix: reformation #10

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056205

 10
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 11
	stiffness matrix reformations = 10
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    2.514453e-04    0.000000e+00 
	   energy              1.238064e+03    1.715156e-05    1.238064e+01 
	   displacement        3.719669e+02    4.111720e-04    3.756603e-04 
Reforming stiffness matrix: reformation #11

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056205

 11
 Nonlinear solution status: time= 0.8
	stiffness updates             = 0
	right hand side evaluations   = 12
	stiffness matrix reformations = 11
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.845164e+04    1.429608e-04    0.000000e+00 
	   energy              1.238064e+03    4.743539e-05    1.238064e+01 
	   displacement        3.719669e+02    9.422993e-05    3.756628e-04 

------- converged at time : 0.8

 ]0;(80%) tempModel.feb - FEBio 2.9.1.0  
AUTO STEPPER: decreasing time step, dt = 0.0953928


===== beginning time step 9 : 0.895393 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1056223

 1
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    1.815348e+03    0.000000e+00 
	   energy              1.126098e+03    6.629003e+01    1.126098e+01 
	   displacement        3.363405e+02    3.363405e+02    3.363405e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1078669

 2
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    1.948173e+00    0.000000e+00 
	   energy              1.126098e+03    5.548588e-01    1.126098e+01 
	   displacement        3.363405e+02    4.532764e+00    3.393411e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1075195

 3
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    2.869969e-01    0.000000e+00 
	   energy              1.126098e+03    5.447873e-02    1.126098e+01 
	   displacement        3.363405e+02    4.608392e-01    3.407785e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1073197

 4
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    1.632977e-01    0.000000e+00 
	   energy              1.126098e+03    4.470216e-02    1.126098e+01 
	   displacement        3.363405e+02    1.039593e+00    3.458765e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1072315

 5
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    6.195628e-02    0.000000e+00 
	   energy              1.126098e+03    1.128813e-02    1.126098e+01 
	   displacement        3.363405e+02    1.317570e+00    3.547719e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071703

 6
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    4.522568e-02    0.000000e+00 
	   energy              1.126098e+03    3.219762e-03    1.126098e+01 
	   displacement        3.363405e+02    7.147419e-02    3.560980e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071343

 7
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 7
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    4.197696e-02    0.000000e+00 
	   energy              1.126098e+03    2.032538e-03    1.126098e+01 
	   displacement        3.363405e+02    2.733883e-02    3.567068e-04 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071307

 8
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 8
	step from line search         = 0.592532
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    4.900220e-02    0.000000e+00 
	   energy              1.126098e+03    2.857672e-03    1.126098e+01 
	   displacement        3.363405e+02    4.253632e-03    3.568416e-04 
Reforming stiffness matrix: reformation #9

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071307

 9
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 11
	stiffness matrix reformations = 9
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    2.378702e-02    0.000000e+00 
	   energy              1.126098e+03    5.657659e-03    1.126098e+01 
	   displacement        3.363405e+02    6.703436e-03    3.569103e-04 
Reforming stiffness matrix: reformation #10

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071307

 10
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 13
	stiffness matrix reformations = 10
	step from line search         = 0.494021
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    2.780122e-02    0.000000e+00 
	   energy              1.126098e+03    4.810793e-03    1.126098e+01 
	   displacement        3.363405e+02    3.456138e-03    3.569484e-04 
Reforming stiffness matrix: reformation #11

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071307

 11
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 15
	stiffness matrix reformations = 11
	step from line search         = 0.622777
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    8.212348e-03    0.000000e+00 
	   energy              1.126098e+03    9.344548e-04    1.126098e+01 
	   displacement        3.363405e+02    3.811569e-04    3.569497e-04 
Reforming stiffness matrix: reformation #12

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071307

 12
 Nonlinear solution status: time= 0.895393
	stiffness updates             = 0
	right hand side evaluations   = 16
	stiffness matrix reformations = 12
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.619656e+04    1.121508e-02    0.000000e+00 
	   energy              1.126098e+03    7.195333e-04    1.126098e+01 
	   displacement        3.363405e+02    2.430382e-04    3.569538e-04 

------- converged at time : 0.895393

 ]0;(90%) tempModel.feb - FEBio 2.9.1.0  
AUTO STEPPER: decreasing time step, dt = 0.0871684


===== beginning time step 10 : 0.982561 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1071325

 1
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    1.158064e+03    0.000000e+00 
	   energy              9.399414e+02    3.831845e+01    9.399414e+00 
	   displacement        2.848087e+02    2.848087e+02    2.848087e-04 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1085311

 2
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    1.502512e+00    0.000000e+00 
	   energy              9.399414e+02    4.332983e-01    9.399414e+00 
	   displacement        2.848087e+02    3.296236e+00    2.907965e-04 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083601

 3
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    3.844409e-01    0.000000e+00 
	   energy              9.399414e+02    7.225414e-02    9.399414e+00 
	   displacement        2.848087e+02    4.725317e-01    2.937906e-04 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083133

 4
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    2.365528e-01    0.000000e+00 
	   energy              9.399414e+02    6.403707e-02    9.399414e+00 
	   displacement        2.848087e+02    1.701375e+00    3.021819e-04 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1082431

 5
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    7.326803e-02    0.000000e+00 
	   energy              9.399414e+02    2.038546e-02    9.399414e+00 
	   displacement        2.848087e+02    1.132854e+00    3.108692e-04 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1082125

 6
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    2.780117e-02    0.000000e+00 
	   energy              9.399414e+02    4.447274e-03    9.399414e+00 
	   displacement        2.848087e+02    1.440637e-01    3.138779e-04 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081909

 7
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 7
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    1.091443e-02    0.000000e+00 
	   energy              9.399414e+02    2.056979e-03    9.399414e+00 
	   displacement        2.848087e+02    2.781008e-02    3.149348e-04 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081909

 8
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 9
	stiffness matrix reformations = 8
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    3.525240e-03    0.000000e+00 
	   energy              9.399414e+02    7.019614e-04    9.399414e+00 
	   displacement        2.848087e+02    8.679476e-03    3.154402e-04 
Reforming stiffness matrix: reformation #9

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081927

 9
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 9
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    7.263913e-04    0.000000e+00 
	   energy              9.399414e+02    2.627521e-04    9.399414e+00 
	   displacement        2.848087e+02    3.694758e-03    3.157154e-04 
Reforming stiffness matrix: reformation #10

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081909

 10
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 11
	stiffness matrix reformations = 10
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    2.229660e-05    0.000000e+00 
	   energy              9.399414e+02    2.085945e-05    9.399414e+00 
	   displacement        2.848087e+02    5.559855e-04    3.158165e-04 
Reforming stiffness matrix: reformation #11

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081909

 11
 Nonlinear solution status: time= 0.982561
	stiffness updates             = 0
	right hand side evaluations   = 12
	stiffness matrix reformations = 11
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            2.230042e+04    4.550160e-07    0.000000e+00 
	   energy              9.399414e+02    4.418224e-08    9.399414e+00 
	   displacement        2.848087e+02    1.026392e-05    3.158267e-04 

------- converged at time : 0.982561

 ]0;(98%) tempModel.feb - FEBio 2.9.1.0  
AUTO STEPPER: decreasing time step, dt = 0.0831584

MUST POINT CONTROLLER: adjusting time step. dt = 0.0174388


===== beginning time step 11 : 1 =====
Reforming stiffness matrix: reformation #1

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1081927

 1
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 2
	stiffness matrix reformations = 1
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    1.001577e+01    0.000000e+00 
	   energy              3.798722e+01    2.439366e-01    3.798722e-01 
	   displacement        1.178387e+01    1.178387e+01    1.178387e-05 
Reforming stiffness matrix: reformation #2

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083493

 2
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 3
	stiffness matrix reformations = 2
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    2.594819e-01    0.000000e+00 
	   energy              3.798722e+01    3.329396e-02    3.798722e-01 
	   displacement        1.178387e+01    2.381143e-02    1.191267e-05 
Reforming stiffness matrix: reformation #3

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083349

 3
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 4
	stiffness matrix reformations = 3
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    2.461196e-02    0.000000e+00 
	   energy              3.798722e+01    5.466944e-03    3.798722e-01 
	   displacement        1.178387e+01    1.159268e-02    1.199214e-05 
Reforming stiffness matrix: reformation #4

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083277

 4
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 5
	stiffness matrix reformations = 4
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    8.357746e-03    0.000000e+00 
	   energy              3.798722e+01    2.882484e-03    3.798722e-01 
	   displacement        1.178387e+01    5.235466e-02    1.228383e-05 
Reforming stiffness matrix: reformation #5

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083331

 5
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 6
	stiffness matrix reformations = 5
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    7.338569e-03    0.000000e+00 
	   energy              3.798722e+01    1.845042e-03    3.798722e-01 
	   displacement        1.178387e+01    8.264195e-02    1.277578e-05 
Reforming stiffness matrix: reformation #6

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083169

 6
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 7
	stiffness matrix reformations = 6
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    5.700700e-04    0.000000e+00 
	   energy              3.798722e+01    1.678084e-04    3.798722e-01 
	   displacement        1.178387e+01    6.990776e-03    1.291127e-05 
Reforming stiffness matrix: reformation #7

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083187

 7
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 8
	stiffness matrix reformations = 7
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    1.516314e-04    0.000000e+00 
	   energy              3.798722e+01    6.121028e-05    3.798722e-01 
	   displacement        1.178387e+01    9.746739e-04    1.295532e-05 
Reforming stiffness matrix: reformation #8

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083187

 8
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 9
	stiffness matrix reformations = 8
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    2.840620e-05    0.000000e+00 
	   energy              3.798722e+01    1.405785e-05    3.798722e-01 
	   displacement        1.178387e+01    1.301412e-04    1.296795e-05 
Reforming stiffness matrix: reformation #9

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083187

 9
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 10
	stiffness matrix reformations = 9
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    1.179752e-06    0.000000e+00 
	   energy              3.798722e+01    1.189886e-06    3.798722e-01 
	   displacement        1.178387e+01    2.085676e-05    1.297147e-05 
Reforming stiffness matrix: reformation #10

===== reforming stiffness matrix:
	Nr of equations ........................... : 25447
	Nr of nonzeroes in stiffness matrix ....... : 1083187

 10
 Nonlinear solution status: time= 1
	stiffness updates             = 0
	right hand side evaluations   = 11
	stiffness matrix reformations = 10
	step from line search         = 1.000000
	convergence norms :     INITIAL         CURRENT         REQUIRED
	   residual            9.260780e+02    9.568330e-11    0.000000e+00 
	   energy              3.798722e+01    7.843945e-10    3.798722e-01 
	   displacement        1.178387e+01    1.063156e-06    1.297212e-05 

------- converged at time : 1

 ]0;(100%) tempModel.feb - FEBio 2.9.1.0  

N O N L I N E A R   I T E R A T I O N   I N F O R M A T I O N

	Number of time steps completed .................... : 11

	Total number of equilibrium iterations ............ : 83

	Average number of equilibrium iterations .......... : 7.54545

	Total number of right hand evaluations ............ : 102

	Total number of stiffness reformations ............ : 83

	Time in linear solver: 0:00:32


 Elapsed time : 0:01:22


 N O R M A L   T E R M I N A T I O N

Waiting for log file...
Proceeding to check log file...10-Aug-2020 10:39:14
------- 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.895393
------- converged at time : 0.982561
------- converged at time : 1
--- Done --- 10-Aug-2020 10:39:14

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

        forceDifference=abs(Fz_final-forceBody);

Use inter/extra-polation to guess displacement

        u=time_mat.*displacementMagnitude;
        f=F_reaction(:,3);
        displacementProposed=interp1(f,u,forceBody,'linear','extrap');
        cFigure; hold on;
        hp(1)=plot(f,u,'r-','LineWidth',lineWidth);
        hp(2)=plot(forceBody,displacementProposed,'b.','MarkerSize',25);
        hl=legend(hp,{'displacement - force FEA','Proposed displacement'},'FontSize',fontSize);
        axis square; axis tight; grid on; box on;
        drawnow;
     displacementDifference=displacementProposed-displacementMagnitude;
     if displacementDifference<-maxProposedDisplacementDifference
         displacementMagnitude=displacementMagnitude-maxProposedDisplacementDifference;
     else
         displacementMagnitude=displacementProposed;
     end
        disp(['Force is: ',num2str(Fz_final),' N. Setting displacement to: ',num2str(displacementMagnitude)]);

        if forceDifference<forceDifferenceTolerance
            break
        end
Force is: -320.272 N. Setting displacement to: -2.6265
    end
end

Visualize force

cFigure; hold on;
hp(1)=plot(time_mat,F_reaction(:,1),'r-','LineWidth',lineWidth);
hp(2)=plot(time_mat,F_reaction(:,2),'g-','LineWidth',lineWidth);
hp(3)=plot(time_mat,F_reaction(:,3),'b-','LineWidth',lineWidth);
hl=legend(hp,{'$F_x$','$F_y$','$F_z$'},'Interpreter','Latex','FontSize',fontSize);
axis square; axis tight; grid on; box on;
drawnow;

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_foot,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_foot,C_energy_foot]=element2patch(E_foot,E_energy(:,:,end));
%     [FE_foot,C_energy_foot]=element2patch(E_foot,E_energy(1:size(E_foot,1),:,1));
indBoundaryFacesFoot=tesBoundary(FE_foot,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_foot(indBoundaryFacesFoot,:),V_def,C_energy_foot(indBoundaryFacesFoot));
hp1=gpatch(FE_foot(indBoundaryFacesFoot,:),V_def,CV,'k',1); %Add graphics object to animate
hp1.FaceColor='Interp';
hp2=gpatch(F_E_sole{1},V_def,'kw','none',0.25); %Add graphics object to animate

axisGeom(gca,fontSize);
colormap(gjet(250)); colorbar;
% caxis([0 max(C_energy_foot)/100]);
caxis([0 maxLevelColorbar_SED]);
axis([min(X_DEF(:)) max(X_DEF(:)) min(Y_DEF(:)) max(Y_DEF(:)) min(Z_DEF(:)) max(Z_DEF(:))]);
view([-40 -40]);
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_foot,C_energy_foot]=element2patch(E_foot,E_energy(:,:,qt));
    CV=faceToVertexMeasure(FE_foot(indBoundaryFacesFoot,:),V_def,C_energy_foot(indBoundaryFacesFoot));

    %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
drawnow;

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