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
- 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 --- 13-Mar-2020 12:22:39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- Writing SMESH file --- 13-Mar-2020 12:22:39 ----> Adding node field ----> Adding facet field ----> Adding holes specification ----> Adding region specification --- Done --- 13-Mar-2020 12:22:39 --- Running TetGen to mesh input boundary--- 13-Mar-2020 12:22:39 Opening /mnt/data/MATLAB/GIBBON/data/temp/temp.smesh. Delaunizing vertices... Delaunay seconds: 0.039103 Creating surface mesh ... Surface mesh seconds: 0.011248 Recovering boundaries... Boundary recovery seconds: 0.027423 Removing exterior tetrahedra ... Spreading region attributes. Exterior tets removal seconds: 0.017467 Recovering Delaunayness... Delaunay recovery seconds: 0.012961 Refining mesh... Refinement seconds: 0.153073 Optimizing mesh... Optimization seconds: 0.009416 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.141619 Total running seconds: 0.412759 Statistics: Input points: 6540 Input facets: 13072 Input segments: 19608 Input holes: 1 Input regions: 1 Mesh points: 10963 Mesh tetrahedra: 49003 Mesh faces: 104542 Mesh faces on exterior boundary: 13072 Mesh faces on input facets: 13072 Mesh edges on input segments: 19608 Steiner points inside domain: 4423 --- Done --- 13-Mar-2020 12:22:40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- Importing TetGen files --- 13-Mar-2020 12:22:40 --- Done --- 13-Mar-2020 12:22:40
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,'bw','b',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,V,'w','none',0.1); 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
%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; %Material section 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 --- 13-Mar-2020 12:22:50
===========================================================================
________ _________ _________ __ _________
| |\ | |\ | |\ | |\ / \\
| ____|| | ____|| | __ || |__|| | ___ ||
| |\___\| | |\___\| | |\_| || \_\| | // \ ||
| || | || | || | || __ | || | ||
| ||__ | ||__ | ||_| || | |\ | || | ||
| |\ | |\ | \\ | || | || | ||
| ___|| | ___|| | ___ || | || | || | ||
| |\__\| | |\__\| | |\__| || | || | || | ||
| || | || | || | || | || | || | ||
| || | ||___ | ||__| || | || | \\__/ ||
| || | |\ | || | || | ||
|___|| |________|| |__________|| |__|| \_________//
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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1005190
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.787056e+04 1.796322e-22 0.000000e+00
energy 1.239875e+03 8.598305e-12 1.239875e+01
displacement 4.124547e+02 4.124547e+02 4.124547e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1005190
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.787056e+04 2.097654e+00 0.000000e+00
energy 1.239875e+03 5.605370e-01 1.239875e+01
displacement 4.124547e+02 4.124547e+02 4.124547e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1006558
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.787056e+04 2.609836e-02 0.000000e+00
energy 1.239875e+03 1.085745e-02 1.239875e+01
displacement 4.124547e+02 4.816696e-02 4.117548e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1006216
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.787056e+04 1.816699e-04 0.000000e+00
energy 1.239875e+03 6.768522e-06 1.239875e+01
displacement 4.124547e+02 2.223353e-03 4.115790e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1006216
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.787056e+04 2.087088e-08 0.000000e+00
energy 1.239875e+03 4.469180e-07 1.239875e+01
displacement 4.124547e+02 7.165760e-05 4.115767e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1006216
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.787348e+04 1.722103e+01 0.000000e+00
energy 1.239561e+03 2.390009e+00 1.239561e+01
displacement 4.096040e+02 4.096040e+02 4.096040e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1009510
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.787348e+04 2.104562e-01 0.000000e+00
energy 1.239561e+03 6.199884e-02 1.239561e+01
displacement 4.096040e+02 1.463780e-01 4.086101e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008916
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.787348e+04 9.370039e-03 0.000000e+00
energy 1.239561e+03 2.760157e-03 1.239561e+01
displacement 4.096040e+02 4.789680e-02 4.072461e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008754
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.787348e+04 1.114058e-03 0.000000e+00
energy 1.239561e+03 2.938874e-04 1.239561e+01
displacement 4.096040e+02 2.838250e-03 4.069964e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
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.787348e+04 1.088413e-03 0.000000e+00
energy 1.239561e+03 3.332791e-05 1.239561e+01
displacement 4.096040e+02 6.424351e-04 4.068862e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
6
Nonlinear solution status: time= 0.3
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.787348e+04 1.459179e-03 0.000000e+00
energy 1.239561e+03 1.019832e-04 1.239561e+01
displacement 4.096040e+02 1.106113e-03 4.068260e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
7
Nonlinear solution status: time= 0.3
stiffness updates = 0
right hand side evaluations = 9
stiffness matrix reformations = 7
step from line search = 0.586716
convergence norms : INITIAL CURRENT REQUIRED
residual 2.787348e+04 8.930748e-04 0.000000e+00
energy 1.239561e+03 5.631585e-05 1.239561e+01
displacement 4.096040e+02 4.927070e-04 4.068242e-04
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
8
Nonlinear solution status: time= 0.3
stiffness updates = 0
right hand side evaluations = 11
stiffness matrix reformations = 8
step from line search = 0.591175
convergence norms : INITIAL CURRENT REQUIRED
residual 2.787348e+04 2.720942e-04 0.000000e+00
energy 1.239561e+03 1.446273e-04 1.239561e+01
displacement 4.096040e+02 5.490122e-04 4.068333e-04
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
9
Nonlinear solution status: time= 0.3
stiffness updates = 0
right hand side evaluations = 13
stiffness matrix reformations = 9
step from line search = 0.500803
convergence norms : INITIAL CURRENT REQUIRED
residual 2.787348e+04 3.000237e-04 0.000000e+00
energy 1.239561e+03 2.713233e-05 1.239561e+01
displacement 4.096040e+02 3.906708e-05 4.068287e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1008574
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.789113e+04 2.250555e+01 0.000000e+00
energy 1.239151e+03 3.103029e+00 1.239151e+01
displacement 4.044517e+02 4.044517e+02 4.044517e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1013092
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.789113e+04 2.083536e-01 0.000000e+00
energy 1.239151e+03 9.373163e-02 1.239151e+01
displacement 4.044517e+02 3.606027e-01 4.018699e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1011922
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.789113e+04 4.392346e-03 0.000000e+00
energy 1.239151e+03 1.888481e-03 1.239151e+01
displacement 4.044517e+02 7.273390e-02 4.004870e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1011490
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.789113e+04 2.718162e-04 0.000000e+00
energy 1.239151e+03 4.739171e-05 1.239151e+01
displacement 4.044517e+02 2.268288e-03 4.002549e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1011490
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.789113e+04 1.870046e-04 0.000000e+00
energy 1.239151e+03 1.374579e-05 1.239151e+01
displacement 4.044517e+02 3.830920e-05 4.002507e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1011490
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.793398e+04 1.365974e+02 0.000000e+00
energy 1.238889e+03 9.492616e+00 1.238889e+01
displacement 3.969408e+02 3.969408e+02 3.969408e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1020904
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.793398e+04 9.601382e-01 0.000000e+00
energy 1.238889e+03 2.350116e-01 1.238889e+01
displacement 3.969408e+02 5.473637e-01 3.944169e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1019230
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.793398e+04 2.801279e-02 0.000000e+00
energy 1.238889e+03 8.768358e-03 1.238889e+01
displacement 3.969408e+02 1.415678e-01 3.924549e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1018222
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.793398e+04 2.549438e-03 0.000000e+00
energy 1.238889e+03 6.214277e-04 1.238889e+01
displacement 3.969408e+02 1.032958e-02 3.919769e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1017862
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.793398e+04 1.342158e-03 0.000000e+00
energy 1.238889e+03 1.274436e-03 1.238889e+01
displacement 3.969408e+02 1.934621e-02 3.916384e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1017862
6
Nonlinear solution status: time= 0.5
stiffness updates = 0
right hand side evaluations = 8
stiffness matrix reformations = 6
step from line search = 0.603135
convergence norms : INITIAL CURRENT REQUIRED
residual 2.793398e+04 2.172798e-04 0.000000e+00
energy 1.238889e+03 3.811654e-04 1.238889e+01
displacement 3.969408e+02 8.634988e-03 3.918352e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1017862
7
Nonlinear solution status: time= 0.5
stiffness updates = 0
right hand side evaluations = 10
stiffness matrix reformations = 7
step from line search = 0.501006
convergence norms : INITIAL CURRENT REQUIRED
residual 2.793398e+04 1.171206e-04 0.000000e+00
energy 1.238889e+03 4.693792e-05 1.238889e+01
displacement 3.969408e+02 2.438081e-04 3.918270e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1017862
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.800767e+04 2.763142e+02 0.000000e+00
energy 1.238453e+03 1.871757e+01 1.238453e+01
displacement 3.880372e+02 3.880372e+02 3.880372e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1037662
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.800767e+04 2.260090e+00 0.000000e+00
energy 1.238453e+03 4.880951e-01 1.238453e+01
displacement 3.880372e+02 1.203140e+00 3.852797e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1034872
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.800767e+04 9.906327e-02 0.000000e+00
energy 1.238453e+03 2.228380e-02 1.238453e+01
displacement 3.880372e+02 2.637873e-01 3.826161e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1032712
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.800767e+04 7.451678e-03 0.000000e+00
energy 1.238453e+03 2.131071e-03 1.238453e+01
displacement 3.880372e+02 4.337167e-02 3.820211e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1031794
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.800767e+04 8.264332e-04 0.000000e+00
energy 1.238453e+03 2.517899e-04 1.238453e+01
displacement 3.880372e+02 6.937007e-03 3.819552e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1031578
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.800767e+04 1.420511e-04 0.000000e+00
energy 1.238453e+03 4.034721e-05 1.238453e+01
displacement 3.880372e+02 2.509630e-03 3.819175e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1031578
7
Nonlinear solution status: time= 0.6
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.800767e+04 3.787369e-07 0.000000e+00
energy 1.238453e+03 4.867610e-07 1.238453e+01
displacement 3.880372e+02 7.388610e-05 3.819164e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1031578
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.812714e+04 5.584429e+02 0.000000e+00
energy 1.237178e+03 2.891242e+01 1.237178e+01
displacement 3.737553e+02 3.737553e+02 3.737553e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1053934
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.812714e+04 3.017690e+00 0.000000e+00
energy 1.237178e+03 7.763288e-01 1.237178e+01
displacement 3.737553e+02 2.153642e+00 3.740582e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1050316
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.812714e+04 1.336491e-01 0.000000e+00
energy 1.237178e+03 3.244602e-02 1.237178e+01
displacement 3.737553e+02 3.676329e-01 3.729038e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1048588
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.812714e+04 1.954225e-02 0.000000e+00
energy 1.237178e+03 7.116544e-03 1.237178e+01
displacement 3.737553e+02 2.175623e-01 3.727095e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047652
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.812714e+04 4.411175e-03 0.000000e+00
energy 1.237178e+03 1.859939e-03 1.237178e+01
displacement 3.737553e+02 7.942120e-02 3.730796e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047364
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.812714e+04 1.848239e-03 0.000000e+00
energy 1.237178e+03 4.077019e-04 1.237178e+01
displacement 3.737553e+02 6.239163e-03 3.732067e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047256
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.812714e+04 9.016213e-04 0.000000e+00
energy 1.237178e+03 1.747055e-04 1.237178e+01
displacement 3.737553e+02 1.548487e-03 3.732262e-04
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047238
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.812714e+04 1.365365e-05 0.000000e+00
energy 1.237178e+03 4.574865e-05 1.237178e+01
displacement 3.737553e+02 1.100289e-03 3.732093e-04
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047238
9
Nonlinear solution status: time= 0.7
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.812714e+04 3.126822e-10 0.000000e+00
energy 1.237178e+03 2.307643e-08 1.237178e+01
displacement 3.737553e+02 4.246053e-06 3.732100e-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 ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1047202
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.831646e+04 7.199969e+02 0.000000e+00
energy 1.236419e+03 3.807546e+01 1.236419e+01
displacement 3.739304e+02 3.739304e+02 3.739304e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1070188
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.831646e+04 2.598768e+00 0.000000e+00
energy 1.236419e+03 6.458493e-01 1.236419e+01
displacement 3.739304e+02 2.484371e+00 3.748844e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1066768
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.831646e+04 1.819910e-01 0.000000e+00
energy 1.236419e+03 3.939338e-02 1.236419e+01
displacement 3.739304e+02 4.209526e-01 3.735262e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1064644
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.831646e+04 6.588946e-02 0.000000e+00
energy 1.236419e+03 1.913890e-02 1.236419e+01
displacement 3.739304e+02 3.241463e-01 3.745058e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1063348
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.831646e+04 3.295739e-02 0.000000e+00
energy 1.236419e+03 4.408418e-03 1.236419e+01
displacement 3.739304e+02 4.416815e-01 3.769015e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062862
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.831646e+04 8.928106e-03 0.000000e+00
energy 1.236419e+03 1.236309e-03 1.236419e+01
displacement 3.739304e+02 5.730830e-02 3.776230e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062646
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.831646e+04 2.330891e-02 0.000000e+00
energy 1.236419e+03 3.325264e-03 1.236419e+01
displacement 3.739304e+02 9.329633e-03 3.778358e-04
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062646
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.831646e+04 1.187611e-02 0.000000e+00
energy 1.236419e+03 2.430642e-03 1.236419e+01
displacement 3.739304e+02 4.247735e-03 3.778612e-04
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062646
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.831646e+04 1.646694e-02 0.000000e+00
energy 1.236419e+03 2.948592e-03 1.236419e+01
displacement 3.739304e+02 3.571670e-03 3.778315e-04
Reforming stiffness matrix: reformation #10
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062646
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.831646e+04 5.754030e-03 0.000000e+00
energy 1.236419e+03 1.690045e-03 1.236419e+01
displacement 3.739304e+02 3.146832e-03 3.778046e-04
Reforming stiffness matrix: reformation #11
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062646
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.831646e+04 6.419186e-03 0.000000e+00
energy 1.236419e+03 1.549591e-03 1.236419e+01
displacement 3.739304e+02 1.796040e-03 3.777765e-04
Reforming stiffness matrix: reformation #12
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062700
12
Nonlinear solution status: time= 0.8
stiffness updates = 0
right hand side evaluations = 13
stiffness matrix reformations = 12
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 2.831646e+04 7.633837e-03 0.000000e+00
energy 1.236419e+03 3.678916e-04 1.236419e+01
displacement 3.739304e+02 1.178137e-03 3.777824e-04
Reforming stiffness matrix: reformation #13
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062700
13
Nonlinear solution status: time= 0.8
stiffness updates = 0
right hand side evaluations = 15
stiffness matrix reformations = 13
step from line search = 0.556838
convergence norms : INITIAL CURRENT REQUIRED
residual 2.831646e+04 1.046965e-02 0.000000e+00
energy 1.236419e+03 9.402932e-04 1.236419e+01
displacement 3.739304e+02 7.069564e-04 3.777679e-04
Reforming stiffness matrix: reformation #14
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062700
14
Nonlinear solution status: time= 0.8
stiffness updates = 0
right hand side evaluations = 16
stiffness matrix reformations = 14
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 2.831646e+04 1.451043e-02 0.000000e+00
energy 1.236419e+03 1.211603e-03 1.236419e+01
displacement 3.739304e+02 6.298333e-04 3.777757e-04
Reforming stiffness matrix: reformation #15
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062700
15
Nonlinear solution status: time= 0.8
stiffness updates = 0
right hand side evaluations = 19
stiffness matrix reformations = 15
step from line search = 0.315941
convergence norms : INITIAL CURRENT REQUIRED
residual 2.831646e+04 7.408438e-03 0.000000e+00
energy 1.236419e+03 8.640738e-04 1.236419e+01
displacement 3.739304e+02 3.148576e-04 3.777645e-04
------- converged at time : 0.8
]0;(80%) tempModel.feb - FEBio 2.9.1.0
AUTO STEPPER: decreasing time step, dt = 0.0818332
===== beginning time step 9 : 0.881833 =====
Reforming stiffness matrix: reformation #1
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1062628
1
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 2
stiffness matrix reformations = 1
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 1.061215e+03 0.000000e+00
energy 8.277608e+02 4.028039e+01 8.277608e+00
displacement 2.492718e+02 2.492718e+02 2.492718e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1082302
2
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 3
stiffness matrix reformations = 2
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 1.682896e+00 0.000000e+00
energy 8.277608e+02 4.443259e-01 8.277608e+00
displacement 2.492718e+02 2.558685e+00 2.516055e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1079098
3
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 4
stiffness matrix reformations = 3
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 2.224807e-01 0.000000e+00
energy 8.277608e+02 4.157959e-02 8.277608e+00
displacement 2.492718e+02 3.341720e-01 2.526942e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1077370
4
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 5
stiffness matrix reformations = 4
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 1.076953e-01 0.000000e+00
energy 8.277608e+02 2.957869e-02 8.277608e+00
displacement 2.492718e+02 5.809898e-01 2.554200e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1076578
5
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 6
stiffness matrix reformations = 5
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 5.069453e-02 0.000000e+00
energy 8.277608e+02 9.503665e-03 8.277608e+00
displacement 2.492718e+02 8.047864e-01 2.608152e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075912
6
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 7
stiffness matrix reformations = 6
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 1.266298e-02 0.000000e+00
energy 8.277608e+02 2.393130e-03 8.277608e+00
displacement 2.492718e+02 1.135219e-01 2.622602e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075678
7
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 8
stiffness matrix reformations = 7
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 5.014789e-03 0.000000e+00
energy 8.277608e+02 8.774460e-05 8.277608e+00
displacement 2.492718e+02 1.727181e-02 2.627329e-04
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075660
8
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 9
stiffness matrix reformations = 8
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 5.720579e-03 0.000000e+00
energy 8.277608e+02 6.785223e-04 8.277608e+00
displacement 2.492718e+02 8.844292e-04 2.628132e-04
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075588
9
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 10
stiffness matrix reformations = 9
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 2.551387e-03 0.000000e+00
energy 8.277608e+02 4.402403e-04 8.277608e+00
displacement 2.492718e+02 8.772928e-04 2.627603e-04
Reforming stiffness matrix: reformation #10
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075588
10
Nonlinear solution status: time= 0.881833
stiffness updates = 0
right hand side evaluations = 12
stiffness matrix reformations = 10
step from line search = 0.559850
convergence norms : INITIAL CURRENT REQUIRED
residual 1.918121e+04 1.257270e-04 0.000000e+00
energy 8.277608e+02 1.580275e-05 8.277608e+00
displacement 2.492718e+02 1.067426e-04 2.627590e-04
------- converged at time : 0.881833
]0;(88%) tempModel.feb - FEBio 2.9.1.0
===== beginning time step 10 : 0.963666 =====
Reforming stiffness matrix: reformation #1
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1075588
1
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 2
stiffness matrix reformations = 1
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 1.205023e+03 0.000000e+00
energy 8.273898e+02 4.099237e+01 8.273898e+00
displacement 2.516370e+02 2.516370e+02 2.516370e-04
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1090222
2
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 3
stiffness matrix reformations = 2
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 1.435996e+00 0.000000e+00
energy 8.273898e+02 4.127033e-01 8.273898e+00
displacement 2.516370e+02 3.455140e+00 2.565281e-04
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1088440
3
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 4
stiffness matrix reformations = 3
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 3.532258e-01 0.000000e+00
energy 8.273898e+02 6.184469e-02 8.273898e+00
displacement 2.516370e+02 3.143129e-01 2.586696e-04
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1087738
4
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 5
stiffness matrix reformations = 4
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 1.760236e-01 0.000000e+00
energy 8.273898e+02 5.391190e-02 8.273898e+00
displacement 2.516370e+02 1.363059e+00 2.653211e-04
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1087216
5
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 6
stiffness matrix reformations = 5
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 5.275117e-02 0.000000e+00
energy 8.273898e+02 1.670607e-02 8.273898e+00
displacement 2.516370e+02 1.102399e+00 2.733665e-04
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086118
6
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 7
stiffness matrix reformations = 6
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 1.788545e-02 0.000000e+00
energy 8.273898e+02 3.172355e-03 8.273898e+00
displacement 2.516370e+02 8.939029e-02 2.754797e-04
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086064
7
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 8
stiffness matrix reformations = 7
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 7.359259e-03 0.000000e+00
energy 8.273898e+02 1.345060e-03 8.273898e+00
displacement 2.516370e+02 1.752152e-02 2.760877e-04
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086082
8
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 9
stiffness matrix reformations = 8
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 2.952413e-03 0.000000e+00
energy 8.273898e+02 4.525459e-04 8.273898e+00
displacement 2.516370e+02 5.663173e-03 2.764006e-04
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086100
9
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 10
stiffness matrix reformations = 9
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 1.882748e-03 0.000000e+00
energy 8.273898e+02 2.730614e-04 8.273898e+00
displacement 2.516370e+02 1.404900e-03 2.765395e-04
Reforming stiffness matrix: reformation #10
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086082
10
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 11
stiffness matrix reformations = 10
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 9.154775e-04 0.000000e+00
energy 8.273898e+02 5.897824e-05 8.273898e+00
displacement 2.516370e+02 5.023820e-04 2.765667e-04
Reforming stiffness matrix: reformation #11
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086082
11
Nonlinear solution status: time= 0.963666
stiffness updates = 0
right hand side evaluations = 12
stiffness matrix reformations = 11
step from line search = 1.000000
convergence norms : INITIAL CURRENT REQUIRED
residual 1.948911e+04 5.248040e-05 0.000000e+00
energy 8.273898e+02 2.875790e-05 8.273898e+00
displacement 2.516370e+02 4.838424e-05 2.765716e-04
------- converged at time : 0.963666
]0;(96%) tempModel.feb - FEBio 2.9.1.0
AUTO STEPPER: decreasing time step, dt = 0.0780714
MUST POINT CONTROLLER: adjusting time step. dt = 0.0363337
===== beginning time step 11 : 1 =====
Reforming stiffness matrix: reformation #1
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1086064
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 3.940237e+03 1.412823e+01 0.000000e+00
energy 1.636941e+02 9.494760e-02 1.636941e+00
displacement 5.103260e+01 5.103260e+01 5.103260e-05
Reforming stiffness matrix: reformation #2
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1090492
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 3.940237e+03 4.038932e-01 0.000000e+00
energy 1.636941e+02 5.308038e-02 1.636941e+00
displacement 5.103260e+01 5.268280e-02 5.139238e-05
Reforming stiffness matrix: reformation #3
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1090078
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 3.940237e+03 6.941697e-02 0.000000e+00
energy 1.636941e+02 1.530452e-02 1.636941e+00
displacement 5.103260e+01 4.497812e-02 5.180234e-05
Reforming stiffness matrix: reformation #4
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089736
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 3.940237e+03 4.368358e-02 0.000000e+00
energy 1.636941e+02 1.012663e-02 1.636941e+00
displacement 5.103260e+01 3.411521e-01 5.351097e-05
Reforming stiffness matrix: reformation #5
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089646
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 3.940237e+03 1.308180e-02 0.000000e+00
energy 1.636941e+02 3.182855e-03 1.636941e+00
displacement 5.103260e+01 1.947570e-01 5.500068e-05
Reforming stiffness matrix: reformation #6
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089430
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 3.940237e+03 3.498226e-03 0.000000e+00
energy 1.636941e+02 8.378021e-04 1.636941e+00
displacement 5.103260e+01 2.906466e-02 5.556552e-05
Reforming stiffness matrix: reformation #7
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089376
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 3.940237e+03 9.317263e-04 0.000000e+00
energy 1.636941e+02 3.014558e-04 1.636941e+00
displacement 5.103260e+01 4.511123e-03 5.576635e-05
Reforming stiffness matrix: reformation #8
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089376
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 3.940237e+03 1.062968e-04 0.000000e+00
energy 1.636941e+02 4.940105e-05 1.636941e+00
displacement 5.103260e+01 1.486396e-03 5.587206e-05
Reforming stiffness matrix: reformation #9
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089376
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 3.940237e+03 7.139726e-07 0.000000e+00
energy 1.636941e+02 5.424882e-07 1.636941e+00
displacement 5.103260e+01 8.153484e-05 5.589392e-05
Reforming stiffness matrix: reformation #10
===== reforming stiffness matrix:
Nr of equations ........................... : 25522
Nr of nonzeroes in stiffness matrix ....... : 1089376
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 3.940237e+03 1.690901e-09 0.000000e+00
energy 1.636941e+02 1.268048e-09 1.636941e+00
displacement 5.103260e+01 1.046078e-07 5.589427e-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 ............ : 88
Average number of equilibrium iterations .......... : 8
Total number of right hand evaluations ............ : 108
Total number of stiffness reformations ............ : 88
Time in linear solver: 0:00:32
Elapsed time : 0:01:21
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...13-Mar-2020 12:24:12
------- 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.881833
------- converged at time : 0.963666
------- converged at time : 1
--- Done --- 13-Mar-2020 12:24:13
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: -321.048 N. Setting displacement to: -2.6247
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/.