DEMO_FEBio_bar_force

Below is a demonstration for: 1) The creation of an FEBio model whereby force is applied to a selection of nodes, in this case to the end of a bar 2) Running an FEBio job with MATLAB 3) Importing FEBio results into MATLAB

Contents

clear all; close all; clc;

Plot settings

figColor='w'; figColorDef='white';
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=1;
edgeColor=0.25*ones(1,3);
edgeWidth=1.5;
markerSize=50;

Control parameters

% path names
filePath=mfilename('fullpath');
savePath=fullfile(fileparts(filePath),'data','temp');

modelName=fullfile(savePath,'tempModel');

%Specifying dimensions and number of elements
sampleWidth=20;
sampleThickness=5;
sampleHeight=5;
pointSpacings=1.*[1 1 1];

numElementsWidth=round(sampleWidth/pointSpacings(1));
numElementsThickness=round(sampleThickness/pointSpacings(2));
numElementsHeight=round(sampleHeight/pointSpacings(3));

forceMagnitude=[0 0 -2e-5];

CREATING MESHED BOX

%Create box 1
boxDim=[sampleWidth sampleThickness sampleHeight]; %Dimensions
boxEl=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
[box1]=hexMeshBox(boxDim,boxEl);
E=box1.E;
V=box1.V;
Fb=box1.Fb;
faceBoundaryMarker=box1.faceBoundaryMarker;

X=V(:,1); Y=V(:,2); Z=V(:,3);
VE=[mean(X(E),2) mean(Y(E),2) mean(Z(E),2)];

elementMaterialIndices=ones(size(E,1),1);
% Plotting boundary surfaces
hf=figuremax(figColor,figColorDef);
title('Model surfaces','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;
patch('Faces',Fb,'Vertices',V,'FaceColor','flat','CData',faceBoundaryMarker,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor);

colormap(jet(6)); colorbar;
set(gca,'FontSize',fontSize);
view(3); axis tight;  axis equal;  grid on;
% camlight headlight;
drawnow;

DEFINE BC's

%Supported nodes
logicRigid=faceBoundaryMarker==1;
Fr=Fb(logicRigid,:);
bcRigidList=unique(Fr(:));

%Prescribed force nodes
logicPrescribe=faceBoundaryMarker==2;
Fr=Fb(logicPrescribe,:);
bcPrescribeList=unique(Fr(:));
bcPrescribeMagnitudes=forceMagnitude(ones(1,numel(bcPrescribeList)),:);

Visualize BC's

hf=figuremax(figColor,figColorDef);
title('Complete model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

patch('Faces',Fb,'Vertices',V,'FaceColor','flat','CData',faceBoundaryMarker,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor);
plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize);
plotV(V(bcPrescribeList,:),'k.','MarkerSize',markerSize);
set(gca,'FontSize',fontSize);

view(3); axis tight;  axis equal;  grid on;
drawnow;

CONSTRUCTING FEB MODEL

FEB_struct.febio_spec.version='1.2';

% Defining file names
FEB_struct.run_filename=[modelName,'.feb']; %FEB file name
FEB_struct.run_logname=[modelName,'.txt']; %FEBio log file name

%Geometry section
FEB_struct.Geometry.Nodes=V;
FEB_struct.Geometry.Elements={E}; %The element sets
FEB_struct.Geometry.ElementType={'hex8'}; %The element types
FEB_struct.Geometry.ElementMat={elementMaterialIndices};

%Material section

%Material 1 uncoupled hyperelastic
c1=1e-3;
k=c1*1000;
Mat1.type='Mooney-Rivlin';
Mat1.props={'c1','c2','k'};
Mat1.vals={c1,0,k};
Mat1.aniso_type='none';

%Collect materials in cell array
FEB_struct.Materials={Mat1};

%Step specific control sections
FEB_struct.Control.AnalysisType='static';
FEB_struct.Control.Properties={'time_steps','step_size',...
    'max_refs','max_ups',...
    'dtol','etol','rtol','lstol'};
FEB_struct.Control.Values={10,0.1,...
    25,0,...
    0.001,0.01,0,0.9};
FEB_struct.Control.TimeStepperProperties={'dtmin','dtmax','max_retries','opt_iter','aggressiveness'};
FEB_struct.Control.TimeStepperValues={1e-5, 0.1, 5, 5, 1};

%Adding BC information
FEB_struct.Boundary.FixList={bcRigidList,bcRigidList,bcRigidList};
FEB_struct.Boundary.FixType={'x','y','z'};

%Loads
FEB_struct.Loads.Force.Nodes={bcPrescribeList,bcPrescribeList,bcPrescribeList};
FEB_struct.Loads.Force.PrescribeType={'x','y','z'};
FEB_struct.Loads.Force.PrescribeValues={bcPrescribeMagnitudes(:,1),bcPrescribeMagnitudes(:,2),bcPrescribeMagnitudes(:,3)};
FEB_struct.Loads.Force.LoadCurveIds=[1,1,1];

%Load curves
FEB_struct.LoadData.LoadCurves.id=[1];
FEB_struct.LoadData.LoadCurves.type={'linear','linear'};
FEB_struct.LoadData.LoadCurves.loadPoints={[0 0;1 1;]};

%Adding output requests
FEB_struct.Output.VarTypes={'displacement','stress','relative volume','shell thickness'};

%Specify log file output
run_node_output_name=[FEB_struct.run_filename(1:end-4),'_node_out.txt'];
FEB_struct.run_output_names={run_node_output_name};
FEB_struct.output_types={'node_data'};
FEB_struct.data_types={'ux;uy;uz'};

SAVING .FEB FILE

FEB_struct.disp_opt=0; %Display waitbars option
febStruct2febFile_v1p2(FEB_struct);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing FEBio XML object --- 23-Sep-2014 12:23:03
Adding Module level
Adding Globals level
Adding Material level
Adding Geometry level
----> Adding node field
----> Adding element field
----> Adding hex8 element entries....
Adding Output level
----> Adding plotfile field
----> Adding logfile field
Adding Boundary level
----> Defining fix type boundary conditions
Adding Loads level
----> Force load
----> Adding nodal forces
----> Adding nodal forces
----> Adding nodal forces
Adding LoadData level
----> Defining load curves
Writing .feb file
--- Done --- 23-Sep-2014 12:23:04

RUNNING FEBIO JOB

% FEBioRunStruct.FEBioPath='C:\Progra~1\FEBio1p8\febio.exe';
% FEBioRunStruct.FEBioPath='C:\Progra~1\FEBio2p0\bin\FEBio2x64.exe';
FEBioRunStruct.run_filename=FEB_struct.run_filename;
FEBioRunStruct.run_logname=FEB_struct.run_logname;
FEBioRunStruct.disp_on=1;
FEBioRunStruct.disp_log_on=1;
% FEBioRunStruct.run_string_quit=run_string_quit;
FEBioRunStruct.t_check=0.25; %Time for checking log file (dont set too small)
FEBioRunStruct.maxtpi=1e99; %Max analysis time
FEBioRunStruct.maxLogCheckTime=3;

%-------------------------------------------------------------------
[rundFlag]=runMonitorFEBio(FEBioRunStruct);%START FEBio NOW!!!!!!!!
%-------------------------------------------------------------------
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- STARTING FEBIO JOB --- 23-Sep-2014 12:23:04
Waiting for log file...
Proceeding to check log file...23-Sep-2014 12:23:04
------- converged at time : 0.1
------- converged at time : 0.170714
------- converged at time : 0.22662
------- converged at time : 0.273871
------- converged at time : 0.313806
------- converged at time : 0.350263
------- converged at time : 0.383545
------- converged at time : 0.413928
------- converged at time : 0.441664
------- converged at time : 0.466984
------- converged at time : 0.490099
------- converged at time : 0.513215
------- converged at time : 0.53633
------- converged at time : 0.559445
------- converged at time : 0.58256
------- converged at time : 0.605675
------- converged at time : 0.628791
------- converged at time : 0.651906
------- converged at time : 0.675021
------- converged at time : 0.698136
------- converged at time : 0.721251
------- converged at time : 0.744367
------- converged at time : 0.767482
------- converged at time : 0.790597
------- converged at time : 0.813712
------- converged at time : 0.836828
------- converged at time : 0.859943
------- converged at time : 0.883058
------- converged at time : 0.906173
------- converged at time : 0.929288
------- converged at time : 0.952404
------- converged at time : 0.975519
------- converged at time : 0.998634
------- converged at time : 1
--- Done --- 23-Sep-2014 12:23:12

IMPORTING NODAL DISPLACEMENT RESULTS

Importing nodal displacements from a log file

[~, N_disp_mat,~]=importFEBio_logfile(FEB_struct.run_output_names{1}); %Nodal displacements

DN=N_disp_mat(:,2:end,end); %Final nodal displacements

CREATING NODE SET IN DEFORMED STATE

V_def=V+DN;
DN_magnitude=sqrt(sum(DN.^2,2));

Plotting the deformed model

[CF]=vertexToFaceMeasure(Fb,DN_magnitude);

hf1=figuremax(figColor,figColorDef);
title('The deformed model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

hps=patch('Faces',Fb,'Vertices',V_def,'FaceColor','flat','CData',CF);

view(3); axis tight;  axis equal;  grid on;
colormap jet; colorbar;
% camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

GIBBON

Kevin M. Moerman (kevinmoerman@hotmail.com)