DEMO_FEBio_cube_multi_generation
Below is a demonstration for: 1) The creation of an FEBio model for spherical indentation 2) The use of multiple steps where contact is switched on and off depending on step 4) Running an FEBio job with MATLAB 5) 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,'cubeTensionMultiGeneration_01'); %Specifying dimensions and number of elements sampleWidth=5; sampleThickness=5; sampleHeight=5; pointSpacings=[1 1 1]; numElementsWidth=round(sampleWidth/pointSpacings(1)); numElementsThickness=round(sampleThickness/pointSpacings(2)); numElementsHeight=round(sampleHeight/pointSpacings(3)); displacementMagnitude=[0 0 2];
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)]; logicMaterial_1=VE(:,1)<0; elementMaterialIndices=logicMaterial_1+1; %Creating faces for plotting [F,CF]=hex2patch(E,elementMaterialIndices);
% 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; % Plotting model hf=figuremax(figColor,figColorDef); title('Model materials','FontSize',fontSize); xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on; patch('Faces',F,'Vertices',V,'FaceColor','flat','CData',CF,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor); colormap jet; colorbar; set(gca,'FontSize',fontSize); view(3); axis tight; axis equal; grid on; % camlight headlight; drawnow;


DEFINE BC's
%Supported nodes logicRigid_X=faceBoundaryMarker==1; Fr=Fb(logicRigid_X,:); bcRigidList_X=unique(Fr(:)); logicRigid_Y=faceBoundaryMarker==3; Fr=Fb(logicRigid_Y,:); bcRigidList_Y=unique(Fr(:)); logicRigid_Z=faceBoundaryMarker==5; Fr=Fb(logicRigid_Z,:); bcRigidList_Z=unique(Fr(:)); %Prescribed displacement nodes logicPrescribe=faceBoundaryMarker==6; Fr=Fb(logicPrescribe,:); bcPrescribeList=unique(Fr(:)); bcPrescribeMagnitudes=displacementMagnitude(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',F,'Vertices',V,'FaceColor','flat','CData',CF,'FaceAlpha',faceAlpha1,'lineWidth',edgeWidth,'edgeColor',edgeColor); plotV(V(bcRigidList_X,:),'r.','MarkerSize',markerSize); plotV(V(bcRigidList_Y,:),'g.','MarkerSize',markerSize); plotV(V(bcRigidList_Z,:),'b.','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='2.0'; FEB_struct.Module.Type='solid'; % Defining file names FEB_struct.run_filename=[modelName,'.feb']; %FEB file name FEB_struct.run_logname=[modelName,'.txt']; %FEBio log file name febMatID=elementMaterialIndices; febMatID(elementMaterialIndices==-2)=1; febMatID(elementMaterialIndices==-3)=2; %Creating FEB_struct 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={febMatID}; FEB_struct.Geometry.ElementsPartName={'Cube'}; % DEFINING MATERIALS %Material 1 deformable block c1=1e-3; k=c1*1e3; FEB_struct.Materials{1}.Type='Mooney-Rivlin'; FEB_struct.Materials{1}.Name='cube_mat'; FEB_struct.Materials{1}.Properties={'c1','c2','k'}; FEB_struct.Materials{1}.Values={c1,0,k}; %Material 2 Homes-Mow compressible multigeneration E=1e-3; E_g=[E/100 E]; v_g=[0.3 0.3]; beta_g=[0.1 0.1]; FEB_struct.Materials{2}.Type='multigeneration'; FEB_struct.Materials{2}.Name='Deformable block'; FEB_struct.Materials{2}.Generation{1}.Solid{1}.Type='Holmes-Mow'; FEB_struct.Materials{2}.Generation{1}.Solid{1}.Properties={'E','v','beta'}; FEB_struct.Materials{2}.Generation{1}.Solid{1}.Values={E_g(1),v_g(1),beta_g(1)}; FEB_struct.Materials{2}.Generation{1}.Properties={'start_time'}; FEB_struct.Materials{2}.Generation{1}.Values={0}; FEB_struct.Materials{2}.Generation{2}.Solid{1}.Type='Holmes-Mow'; FEB_struct.Materials{2}.Generation{2}.Solid{1}.Properties={'E','v','beta'}; FEB_struct.Materials{2}.Generation{2}.Solid{1}.Values={E_g(2),v_g(2),beta_g(2)}; FEB_struct.Materials{2}.Generation{2}.Properties={'start_time'}; FEB_struct.Materials{2}.Generation{2}.Values={1}; %Step specific control sections FEB_struct.Step{1}.Control.AnalysisType='static'; FEB_struct.Step{1}.Control.Properties={'time_steps','step_size',... 'max_refs','max_ups',... 'dtol','etol','rtol','lstol'}; FEB_struct.Step{1}.Control.Values={20,0.05,... 25,0,... 0.001,0.01,0,0.9}; FEB_struct.Step{1}.Control.TimeStepperProperties={'dtmin','dtmax','max_retries','opt_iter','aggressiveness'}; FEB_struct.Step{1}.Control.TimeStepperValues={1e-4,0.05,5,10,1}; FEB_struct.Step{2}.Control=FEB_struct.Step{1}.Control; %Defining node sets FEB_struct.Geometry.NodeSet{1}.Set=bcRigidList_X; FEB_struct.Geometry.NodeSet{1}.Name='bcRigidList_X'; FEB_struct.Geometry.NodeSet{2}.Set=bcRigidList_Y; FEB_struct.Geometry.NodeSet{2}.Name='bcRigidList_Y'; FEB_struct.Geometry.NodeSet{3}.Set=bcRigidList_Z; FEB_struct.Geometry.NodeSet{3}.Name='bcRigidList_Z'; % FEB_struct.Geometry.NodeSet{4}.Set=bcPrescribeList; % FEB_struct.Geometry.NodeSet{4}.Name='bcPrescribeList'; %Adding BC information FEB_struct.Boundary.Fix{1}.bc='x'; FEB_struct.Boundary.Fix{1}.SetName=FEB_struct.Geometry.NodeSet{1}.Name; FEB_struct.Boundary.Fix{2}.bc='y'; FEB_struct.Boundary.Fix{2}.SetName=FEB_struct.Geometry.NodeSet{2}.Name; FEB_struct.Boundary.Fix{3}.bc='z'; FEB_struct.Boundary.Fix{3}.SetName=FEB_struct.Geometry.NodeSet{3}.Name; %Step specific BC's FEB_struct.Step{1}.Boundary.Prescribe{1}.Set=bcPrescribeList; FEB_struct.Step{1}.Boundary.Prescribe{1}.bc='z'; FEB_struct.Step{1}.Boundary.Prescribe{1}.lc=1; FEB_struct.Step{1}.Boundary.Prescribe{1}.nodeScale=bcPrescribeMagnitudes(:,3); FEB_struct.Step{2}.Boundary.Prescribe{1}.Set=bcPrescribeList; FEB_struct.Step{2}.Boundary.Prescribe{1}.bc='z'; FEB_struct.Step{2}.Boundary.Prescribe{1}.lc=1; FEB_struct.Step{2}.Boundary.Prescribe{1}.nodeScale=bcPrescribeMagnitudes(:,3); %Load curves FEB_struct.LoadData.LoadCurves.id=1; FEB_struct.LoadData.LoadCurves.type={'linear'}; FEB_struct.LoadData.LoadCurves.loadPoints={[0 0;1 1; 2 0;];}; %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(FEB_struct);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- Writing FEBio XML object --- 18-Aug-2014 21:28:11 Adding Module level Adding Globals level Adding Material level Adding Geometry level ----> Adding node field ----> Adding element field ----> Adding hex8 element entries.... ----> Adding NodeSet field Adding Boundary level ----> Defining fix type boundary conditions Adding LoadData level ----> Defining load curves Adding Step level ----> Adding Module field ----> Adding Boundary field ----> Defining prescribe type boundary conditions ----> Adding Control field Adding Step level ----> Adding Module field ----> Adding Boundary field ----> Defining prescribe type boundary conditions ----> Adding Control field Adding Output level ----> Adding plotfile field ----> Adding logfile field Writing .feb file --- Done --- 18-Aug-2014 21:28:11
RUNNING FEBIO JOB
FEBioRunStruct.FEBioPath='C:\Program Files\febio-2.1.0\bin\FEBio2.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.runMode='external';%'internal'; FEBioRunStruct.t_check=0.25; %Time for checking log file (dont set too small) FEBioRunStruct.maxtpi=1e99; %Max analysis time FEBioRunStruct.maxLogCheckTime=3; %Max log file checking time [runFlag]=runMonitorFEBio(FEBioRunStruct);%START FEBio NOW!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% --- STARTING FEBIO JOB --- 18-Aug-2014 21:28:11 Waiting for log file... Proceeding to check log file...18-Aug-2014 21:28:11 ------- converged at time : 0.05 ------- converged at time : 0.1 ------- converged at time : 0.15 ------- converged at time : 0.2 ------- converged at time : 0.25 ------- converged at time : 0.3 ------- converged at time : 0.35 ------- converged at time : 0.4 ------- converged at time : 0.45 ------- converged at time : 0.5 ------- converged at time : 0.55 ------- converged at time : 0.6 ------- converged at time : 0.65 ------- converged at time : 0.7 ------- converged at time : 0.75 ------- converged at time : 0.8 ------- converged at time : 0.85 ------- converged at time : 0.9 ------- converged at time : 0.95 ------- converged at time : 1 ------- converged at time : 1.05 ------- converged at time : 1.1 ------- converged at time : 1.15 ------- converged at time : 1.2 ------- converged at time : 1.25 ------- converged at time : 1.3 ------- converged at time : 1.35 ------- converged at time : 1.4 ------- converged at time : 1.45 ------- converged at time : 1.5 ------- converged at time : 1.55 ------- converged at time : 1.6 ------- converged at time : 1.65 ------- converged at time : 1.7 ------- converged at time : 1.75 ------- converged at time : 1.8 ------- converged at time : 1.85 ------- converged at time : 1.9 ------- converged at time : 1.95 ------- converged at time : 2 --- Done --- 18-Aug-2014 21:28:13
if runFlag==1 %i.e. a succesful run % 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; end

GIBBON
Kevin M. Moerman (kevinmoerman@hotmail.com)