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)