DEMO_febio_0050_foot_insole_01
Below is a demonstration for:
- Building triangulated surface geometry for a foot
- Meshing the foot using tetrahedral elements
- Building surface model of an insole and meshing it with quads
- Extruding the surface to a thickened layer of hexahedral elements
- Defining the boundary conditions
- Coding the febio structure
- Running the model
- Importing and visualizing results
Contents
- Keywords
- Control parameters
- Import surface models
- Subdivide skin mesh and smooth
- Compute mesh derived parameters
- Reorient surfaces
- Cut skin surface
- Scale foot to desired size
- Close over top of skin
- Joining surface features
- Merge shared nodes
- Mesh foot with tetrahedral elements
- Build insole
- Joining node sets
- Define contact surfaces
- Define boundary conditions
- Set-up materials
- Convert list of materials to element id
- Defining the FEBio input structure
- Running the FEBio analysis
- Export input file
- Run febio
- Set displacement magnitude in input file structure
- Export input file
- Run febio
- Importing rigid body reaction forces from a log file
- Use inter/extra-polation to guess displacement
- Visualize force
- Plot animation
Keywords
- febio_spec version 2.5
- febio, FEBio
- shoe insole
- contact, sliding, friction
- tetrahedral elements, tet4
- hexahedral elements, hex8
- static, solid
- hyperelastic, Ogden
- displacement logfile
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/.