DEMO_FEBio_iFEA_uniaxial_01

Below is a demonstration for: 1) Inverse FEA based material parameter optimisation

Contents

clear; close all; clc;
warning off;

Plot settings

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

Control parameters

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

modelName=fullfile(savePath,'iFEA_tempModel');

%Specifying dimensions and number of elements
sampleWidth=10;
sampleThickness=10;
sampleHeight=10;
pointSpacings=2*ones(1,3);
initialArea=sampleWidth*sampleThickness;

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

stretchLoad=0.7;
displacementMagnitude=[0 0 (stretchLoad*sampleHeight)-sampleHeight];

%True material parameter set
k_factor=1e4;
c1_true=1e-3;
m1_true=12;
k_true=c1_true*k_factor;

%Initial material parameter set
c1_ini=c1_true./2;
m1_ini=m1_true+7;
k_ini=c1_ini*k_factor;
P=[c1_ini m1_ini];

SIMULATE EXPERIMENTAL DATA

%Basic set
stress_cauchy_exp=[0;-0.0422226256000000;-0.0811346871800000;-0.119872916800000;-0.161466624000000;-0.209098742000000;-0.266409832800000;-0.337879334400000;-0.429344276800000;-0.548728823600000;-0.707119980000000];
stretch_exp=[1;0.970000000000000;0.940000000000000;0.910000000000000;0.880000000000000;0.850000000000000;0.820000000000000;0.790000000000000;0.760000000000000;0.730000000000000;0.700000000000000];

%Interpolate to higher sampling
n=100;
stretch_exp_n=linspace(1,stretchLoad,n);
stress_cauchy_exp_n = interp1(stretch_exp,stress_cauchy_exp,stretch_exp_n,'cubic');

%Override variables
stress_cauchy_exp=stress_cauchy_exp_n;
stretch_exp=stretch_exp_n;

%Add noise
stdNoise=0.01; %Standard deviation in units of stress
stress_cauchy_exp_n=stress_cauchy_exp_n+stdNoise.*randn(size(stress_cauchy_exp_n));

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',faceAlpha2,'lineWidth',edgeWidth,'edgeColor',edgeColor);

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

DEFINE BC's

%Define supported node sets
logicFace=faceBoundaryMarker==1;
Fr=Fb(logicFace,:);
bcSupportList_X=unique(Fr(:));

logicFace=faceBoundaryMarker==3;
Fr=Fb(logicFace,:);
bcSupportList_Y=unique(Fr(:));

logicFace=faceBoundaryMarker==5;
Fr=Fb(logicFace,:);
bcSupportList_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',Fb,'Vertices',V,'FaceColor','flat','CData',faceBoundaryMarker,'FaceAlpha',faceAlpha2,'lineWidth',edgeWidth,'edgeColor',edgeColor);
plotV(V(bcSupportList_X,:),'r.','MarkerSize',markerSize);
plotV(V(bcSupportList_Y,:),'g.','MarkerSize',markerSize);
plotV(V(bcSupportList_Z,:),'b.','MarkerSize',markerSize);
plotV(V(bcPrescribeList,:),'k.','MarkerSize',markerSize);
set(gca,'FontSize',fontSize);

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

CONSTRUCTING FEB MODEL

FEB_struct.febio_spec.version='2.0';

% 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};
FEB_struct.Geometry.ElementsPartName={'Block'};

%Material section

%Material 1 uncoupled hyperelastic
FEB_struct.Materials{1}.Type='Ogden';
FEB_struct.Materials{1}.Name='Block_mat';
FEB_struct.Materials{1}.Properties={'c1','m1','k'};
FEB_struct.Materials{1}.Values={c1_ini,m1_ini,k_ini};

%Control section
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={20,0.05,...
    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-4,0.05,5,10,1};

%Defining node sets
FEB_struct.Geometry.NodeSet{1}.Set=bcSupportList_X;
FEB_struct.Geometry.NodeSet{1}.Name='bcSupportList_X';
FEB_struct.Geometry.NodeSet{2}.Set=bcSupportList_Y;
FEB_struct.Geometry.NodeSet{2}.Name='bcSupportList_Y';
FEB_struct.Geometry.NodeSet{3}.Set=bcSupportList_Z;
FEB_struct.Geometry.NodeSet{3}.Name='bcSupportList_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;

%Prescribed BC's
FEB_struct.Boundary.Prescribe{1}.Set=bcPrescribeList;
FEB_struct.Boundary.Prescribe{1}.bc='z';
FEB_struct.Boundary.Prescribe{1}.lc=1;
FEB_struct.Boundary.Prescribe{1}.nodeScale=displacementMagnitude(ones(numel(bcPrescribeList),1),3);
FEB_struct.Boundary.Prescribe{1}.Type='relative';

%Load curves
FEB_struct.LoadData.LoadCurves.id=1;
FEB_struct.LoadData.LoadCurves.type={'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_disp_output_name=[FEB_struct.run_filename(1:end-4),'_node_out.txt'];
run_force_output_name=[FEB_struct.run_filename(1:end-4),'_force_out.txt'];
FEB_struct.run_output_names={run_disp_output_name,run_force_output_name};
FEB_struct.output_types={'node_data','node_data'};
FEB_struct.data_types={'ux;uy;uz','Rx;Ry;Rz'};

SAVING .FEB FILE

FEB_struct.disp_opt=0; %Display waitbars option
febStruct2febFile(FEB_struct);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing FEBio XML object --- 23-Sep-2014 12:38:31
Adding Module level
Adding Control 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
----> Defining prescribe type boundary conditions
Adding LoadData level
----> Defining load curves
Adding Output level
----> Adding plotfile field
----> Adding logfile field
Writing .feb file
--- Done --- 23-Sep-2014 12:38:31

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 --- 23-Sep-2014 12:38:31
Waiting for log file...
Proceeding to check log file...23-Sep-2014 12:38:32
------- converged at time : 0.05
------- converged at time : 0.1
------- converged at time : 0.125
------- converged at time : 0.155
------- converged at time : 0.189
------- converged at time : 0.221422
------- converged at time : 0.239391
------- converged at time : 0.253618
------- converged at time : 0.274999
------- converged at time : 0.302104
------- converged at time : 0.333788
------- converged at time : 0.342625
------- converged at time : 0.348251
------- converged at time : 0.362751
------- converged at time : 0.384352
------- converged at time : 0.411632
------- converged at time : 0.443456
------- converged at time : 0.478916
------- converged at time : 0.515162
------- converged at time : 0.554158
------- converged at time : 0.591345
------- converged at time : 0.631094
------- converged at time : 0.672054
------- converged at time : 0.693437
------- converged at time : 0.712967
------- converged at time : 0.73859
------- converged at time : 0.746215
------- converged at time : 0.753839
------- converged at time : 0.769939
------- converged at time : 0.792819
------- converged at time : 0.821123
------- converged at time : 0.853766
------- converged at time : 0.88988
------- converged at time : 0.928772
------- converged at time : 0.969885
------- converged at time : 1
--- Done --- 23-Sep-2014 12:38:35

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;

IMPORTING NODAL FORCES

Importing nodal forces from a log file

[time_mat, N_force_mat,~]=importFEBio_logfile(FEB_struct.run_output_names{2}); %Nodal displacements

FZ_set=N_force_mat(bcPrescribeList,end,:); %Final nodal displacements
%Get Z forces
FZ=sum(N_force_mat(bcPrescribeList,end,:),1);
FZ=[0; FZ(:)]; %Mean top surface nodal forces

%Derive applied stretch
DZ_set=N_disp_mat(bcPrescribeList,end,:); %Final nodal displacements
DZ_set=mean(DZ_set,1);
stretch_sim=(DZ_set+sampleHeight)./sampleHeight;
stretch_sim=[1; stretch_sim(:)];

%Derive simulated Cauchy stress (alternatively import stress and take the mean)
currentArea=initialArea./stretch_sim;
stress_cauchy_sim=FZ./currentArea; %Cauchy stress
stress_cauchy_sim=stress_cauchy_sim.*1e3; %Scale to kPa

%Interpolate experiment onto simulated points
stress_cauchy_exp_sim = interp1(stretch_exp,stress_cauchy_exp,stretch_sim,'cubic');
hf1=figuremax(figColor,figColorDef);
title('Stretch stress curves, initial differences between simulation and experiment','FontSize',fontSize);
xlabel('\lambda Stretch [.]','FontSize',fontSize); ylabel('\sigma Cauchy stress [kPa]','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

H(1)=plot(stretch_sim,stress_cauchy_sim,'r.-','lineWidth',lineWidth,'markerSize',markerSize);
H(2)=plot(stretch_exp,stress_cauchy_exp,'g-','lineWidth',lineWidth);
H(3)=plot(stretch_sim,stress_cauchy_exp_sim,'k.','markerSize',markerSize);

legend(H,{'Simulation','Experiment','Experiment@simulation'},'Location','northwest');
view(2); axis tight;  grid on;
set(gca,'FontSize',fontSize);
drawnow;

CREATE STRUCTURES FOR OPTIMISATION

mat_struct.id=1;
mat_struct.par_names={'c1','m1','k'};
mat_struct.par_values={c1_ini m1_ini k_ini};

% docNode=set_mat_par_FEBIO(FEB_struct.run_filename,FEB_struct.run_filename,{mat_struct});

FEBioRunStruct.disp_on=0;
FEBioRunStruct.disp_log_on=0;

%What should be known to the objective function:
objectiveStruct.bcPrescribeList=bcPrescribeList;
objectiveStruct.stretch_exp=stretch_exp;
objectiveStruct.stress_cauchy_exp=stress_cauchy_exp;
objectiveStruct.FEBioRunStruct=FEBioRunStruct;
objectiveStruct.FEB_struct=FEB_struct;
objectiveStruct.mat_struct=mat_struct;
objectiveStruct.k_factor=k_factor;
objectiveStruct.P_ini=P;
objectiveStruct.initialArea=initialArea;
objectiveStruct.sampleHeight=sampleHeight;

%Optimisation settings
maxNumberIterations=50; %Maximum number of iterations
maxNumberFunctionEvaluations=maxNumberIterations; %Maximum number of function evaluations, N.B. multiple evaluations are used per iteration
functionTolerance=1e-3; %Tolerance on objective function value
parameterTolerance=1e-3; %Tolerance on parameter variation
displayTypeIterations='iter';
OPT_options = optimset('MaxFunEvals',maxNumberFunctionEvaluations,...
                       'MaxIter',maxNumberIterations,...
                       'TolFun',functionTolerance,...
                       'TolX',parameterTolerance,...
                       'Display',displayTypeIterations);
objectiveStruct.method=1;

%File names of output files
objectiveStruct.run_output_names=FEB_struct.run_output_names;

PARAMETER BOUNDS

Ps.f=1;
Ps.t=0.9;
Ps.ub=[0.1  22];
Ps.lb=[0  2];
Ps.c=P;
objectiveStruct.Ps=Ps;

Pb=zeros(size(P));
for q=1:1:numel(P);
    Psn=objectiveStruct.Ps;
    Psn.ub=objectiveStruct.Ps.ub(q);
    Psn.lb=objectiveStruct.Ps.lb(q);
    Psn.c=objectiveStruct.Ps.c(q);
    Pb(q)=inv_parbound(P(q),Psn);
end

STARTING OPTIMISATION

switch objectiveStruct.method
    case 1 %fminsearch and Nelder-Mead
        [Pb_opt,OPT_out.fval,OPT_out.exitflag,OPT_out.output]= fminsearch(@(Pb) obj_DEMO_FEBio_iFEA_uniaxial_01(Pb,objectiveStruct),Pb,OPT_options);
    case 2 %lsqnonlin and Levenberg-Marquardt
        OPT_options = optimset(OPT_options,'Algorithm',{'levenberg-marquardt',.01}); %Specifically setting algorithm
        [Pb_opt,OPT_out.resnorm,OPT_out.residual]= lsqnonlin(@(Pb) obj_DEMO_FEBio_iFEA_uniaxial_01(Pb,objectiveStruct),Pb,[],[],OPT_options);
end

[Fopt,OPT_stats_out]=obj_DEMO_FEBio_iFEA_uniaxial_01(Pb_opt,objectiveStruct);

type(fullfile(fileparts(filePath),'obj_DEMO_FEBio_iFEA_uniaxial_01'))
SETTING MATERIAL PARAMETERS...
P=4.9999999999999491e-04,1.9000000000000000e+01
 
 Iteration   Func-count     min f(x)         Procedure
     0            1         0.429713         
SETTING MATERIAL PARAMETERS...
P=4.7620536132347293e-04,1.9000000000000000e+01
SETTING MATERIAL PARAMETERS...
P=4.9999999999999491e-04,1.9449093797277914e+01
     1            3         0.413514         initial simplex
SETTING MATERIAL PARAMETERS...
P=5.2629737880616201e-04,1.9449093797277914e+01
SETTING MATERIAL PARAMETERS...
P=5.5551423447968156e-04,1.9628952188266517e+01
     2            5         0.393962         reflect
SETTING MATERIAL PARAMETERS...
P=5.2629737880616201e-04,1.9786088930837408e+01
SETTING MATERIAL PARAMETERS...
P=5.0632487511882989e-04,1.9241562013762504e+01
     3            7         0.393962         contract inside
SETTING MATERIAL PARAMETERS...
P=5.3330966145687662e-04,1.9241562013762504e+01
SETTING MATERIAL PARAMETERS...
P=5.5168596762143480e-04,1.9125498549088075e+01
     4            9         0.366009         expand
SETTING MATERIAL PARAMETERS...
P=5.7548107875163290e-04,1.9349146721670991e+01
SETTING MATERIAL PARAMETERS...
P=5.2200757300061018e-04,1.9269218372207050e+01
     5           11         0.366009         contract inside
SETTING MATERIAL PARAMETERS...
P=5.4697419272429777e-04,1.8899034901191282e+01
SETTING MATERIAL PARAMETERS...
P=5.5793399625689852e-04,1.8534102635957460e+01
     6           13         0.319506         expand
SETTING MATERIAL PARAMETERS...
P=5.9197079666493133e-04,1.8309394330674259e+01
SETTING MATERIAL PARAMETERS...
P=6.3448928245363190e-04,1.7570992694693725e+01
     7           15         0.229257         expand
SETTING MATERIAL PARAMETERS...
P=6.4276740967252890e-04,1.6245406646358191e+01
     8           16         0.229257         reflect
SETTING MATERIAL PARAMETERS...
P=7.4653288508632623e-04,1.4047334211553881e+01
SETTING MATERIAL PARAMETERS...
P=8.9835624098952857e-04,9.9933013770713490e+00
     9           18         0.224649         reflect
SETTING MATERIAL PARAMETERS...
P=7.3538992925490938e-04,1.6070631271161130e+01
SETTING MATERIAL PARAMETERS...
P=7.9248633887606947e-04,1.5980465618454014e+01
    10           20          0.12547         reflect
SETTING MATERIAL PARAMETERS...
P=8.9024057923483559e-04,1.1234375125658248e+01
SETTING MATERIAL PARAMETERS...
P=6.8358720625516511e-04,1.6550951256273557e+01
    11           22          0.12547         contract inside
SETTING MATERIAL PARAMETERS...
P=6.7423207871725741e-04,1.7774272576293967e+01
SETTING MATERIAL PARAMETERS...
P=7.2704217290290869e-04,1.5295687489187756e+01
    12           24          0.12547         contract inside
SETTING MATERIAL PARAMETERS...
P=7.8592254118038944e-04,1.4666019801156949e+01
SETTING MATERIAL PARAMETERS...
P=8.4950684874780531e-04,1.3461256753090328e+01
    13           26         0.054901         expand
SETTING MATERIAL PARAMETERS...
P=8.6092502104772044e-04,1.4465022947340248e+01
SETTING MATERIAL PARAMETERS...
P=9.4822619506012011e-04,1.4011630459308730e+01
    14           28        0.0534792         reflect
SETTING MATERIAL PARAMETERS...
P=1.0215679714094983e-03,1.1351692332346587e+01
SETTING MATERIAL PARAMETERS...
P=1.2683051548957602e-03,8.9217804652593884e+00
    15           30       0.00627733         reflect
SETTING MATERIAL PARAMETERS...
P=1.0381231353120358e-03,1.2457262022400922e+01
    16           31       0.00627733         reflect
SETTING MATERIAL PARAMETERS...
P=1.2809838379581884e-03,9.3716800795088044e+00
SETTING MATERIAL PARAMETERS...
P=9.3781425027776317e-04,1.3234213998883277e+01
    17           33       0.00627733         contract inside
SETTING MATERIAL PARAMETERS...
P=9.2428211280556893e-04,1.2141281624299905e+01
SETTING MATERIAL PARAMETERS...
P=9.5033638079548572e-04,1.2220403015174242e+01
    18           35       0.00627733         contract outside
SETTING MATERIAL PARAMETERS...
P=1.0364435206743018e-03,1.0357464532382828e+01
SETTING MATERIAL PARAMETERS...
P=9.6066941233802217e-04,1.2516304753461913e+01
    19           37        0.0024727         contract inside
SETTING MATERIAL PARAMETERS...
P=1.0335172315461161e-03,1.1646753475466559e+01
SETTING MATERIAL PARAMETERS...
P=1.0808158485901165e-03,1.1361489855456643e+01
    20           39        0.0022504         reflect
SETTING MATERIAL PARAMETERS...
P=9.7122954606024270e-04,1.2809989074825380e+01
SETTING MATERIAL PARAMETERS...
P=1.0085005930330516e-03,1.1715867309363617e+01
    21           41       0.00147326         contract inside
SETTING MATERIAL PARAMETERS...
P=1.0890843426852093e-03,1.0857069116590949e+01
SETTING MATERIAL PARAMETERS...
P=9.8984882241313628e-04,1.2099229068007528e+01
SETTING MATERIAL PARAMETERS...
P=1.0208557261586538e-03,1.1681300994118347e+01
SETTING MATERIAL PARAMETERS...
P=9.8400427755340205e-04,1.2116546084534907e+01
    22           45       0.00080353         shrink
SETTING MATERIAL PARAMETERS...
P=9.9576318728295541e-04,1.2081910582805021e+01
SETTING MATERIAL PARAMETERS...
P=9.8951434410832770e-04,1.2264880448801430e+01
    23           47      4.84663e-05         reflect
SETTING MATERIAL PARAMETERS...
P=1.0335172315461161e-03,1.1646753475466559e+01
SETTING MATERIAL PARAMETERS...
P=9.9593251623556487e-04,1.1999016952617250e+01
SETTING MATERIAL PARAMETERS...
P=1.0081533972502777e-03,1.1881488490712925e+01
SETTING MATERIAL PARAMETERS...
P=9.8984882241313628e-04,1.2099229068007528e+01
    24           51      4.84663e-05         shrink
 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.000048 

SETTING MATERIAL PARAMETERS...
P=9.9576318728295541e-04,1.2081910582805021e+01

function [Fopt,OPT_stats_out]=obj_DEMO_FEBio_iFEA_uniaxial_01(Pb,objectiveStruct)

%% PARAMETER BOUNDS

P=zeros(size(Pb));
for q=1:1:numel(P);
    Psn=objectiveStruct.Ps;
    Psn.ub=objectiveStruct.Ps.ub(q);
    Psn.lb=objectiveStruct.Ps.lb(q);
    Psn.c=objectiveStruct.Ps.c(q);
    P(q)=parbound(Pb(q),Psn);
end

%% SETTING MATERIAL PARAMETERS

%Acces material parameters
mat_struct=objectiveStruct.mat_struct;
mat_struct.par_values={P(1) P(2) P(1)*objectiveStruct.k_factor}; 

disp('SETTING MATERIAL PARAMETERS...');
disp_text=sprintf('%6.16e,',P); disp_text=disp_text(1:end-1);
disp(['P=',disp_text]);

%Assign material parameters
set_mat_par_FEBIO(objectiveStruct.FEB_struct.run_filename,objectiveStruct.FEB_struct.run_filename,{mat_struct});

%% START FEBio NOW

[runFlag]=runMonitorFEBio(objectiveStruct.FEBioRunStruct);

bcPrescribeList=objectiveStruct.bcPrescribeList;
sampleHeight=objectiveStruct.sampleHeight;
initialArea=objectiveStruct.initialArea;
stretch_exp=objectiveStruct.stretch_exp;
stress_cauchy_exp=objectiveStruct.stress_cauchy_exp;

if runFlag==1  
  
    
    %Importing displacement
    [~,N_disp_mat,~]=importFEBio_logfile(objectiveStruct.FEB_struct.run_output_names{1}); %Nodal displacements

    %Importing force
    [~,N_force_mat,~]=importFEBio_logfile(objectiveStruct.FEB_struct.run_output_names{2}); %Nodal displacements
    
    %Get Z forces
    FZ=sum(N_force_mat(bcPrescribeList,end,:),1);
    FZ=[0; FZ(:)]; %Mean top surface nodal forces
    
    %Derive applied stretch
    DZ_set=N_disp_mat(bcPrescribeList,end,:); %Final nodal displacements
    DZ_set=mean(DZ_set,1);
    stretch_sim=(DZ_set+sampleHeight)./sampleHeight;
    stretch_sim=[1; stretch_sim(:)];
    
    %Derive simulated Cauchy stress 
    currentArea=initialArea./stretch_sim;
    stress_cauchy_sim=FZ./currentArea; %Cauchy stress
    stress_cauchy_sim=stress_cauchy_sim.*1e3; %Scale to kPa
    
    %Interpolate experiment onto simulated points
    stress_cauchy_sim_exp = interp1(stretch_sim,stress_cauchy_sim,stretch_exp,'cubic');
    
    %Derive Fopt
    stressDev=stress_cauchy_exp-stress_cauchy_sim_exp;  
       
    switch objectiveStruct.method
        case 1
            Fopt=sum((stressDev).^2); %Sum of squared differences
        case 2
            Fopt=(stressDev).^2; %Squared differences
    end
    
    OPT_stats_out.stress_cauchy_sim=stress_cauchy_sim;
    OPT_stats_out.stretch_sim=stretch_sim;
    OPT_stats_out.stressDev=stressDev;
    OPT_stats_out.Fopt=Fopt;
    
else %Output NaN
    switch objectiveStruct.method
        case 1
            Fopt=NaN; 
        case 2
            Fopt=NaN(size(stress_cauchy_exp)); %Squared differences
    end
    OPT_stats_out=[];
end

%%


P_opt=zeros(size(P));
for q=1:1:numel(P);
    Psn=objectiveStruct.Ps;
    Psn.ub=objectiveStruct.Ps.ub(q);
    Psn.lb=objectiveStruct.Ps.lb(q);
    Psn.c=objectiveStruct.Ps.c(q);
    P_opt(q)=parbound(Pb_opt(q),Psn);
end

disp_text=sprintf('%6.16e,',P_opt); disp_text=disp_text(1:end-1);
disp(['P_opt=',disp_text]);
P_opt=9.9576318728295541e-04,1.2081910582805021e+01
hf1=figuremax(figColor,figColorDef);
title('Stretch stress curves optimised','FontSize',fontSize);
xlabel('\lambda Stretch [.]','FontSize',fontSize); ylabel('\sigma Cauchy stress [kPa]','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

H(1)=plot(OPT_stats_out.stretch_sim,OPT_stats_out.stress_cauchy_sim,'r.-','lineWidth',lineWidth,'markerSize',markerSize);
H(2)=plot(stretch_exp,stress_cauchy_exp,'g-','lineWidth',lineWidth);

legend(H,{'Simulation','Experiment'},'Location','northwest');
view(2); axis tight;  grid on;
set(gca,'FontSize',fontSize);
drawnow;

GIBBON

Kevin M. Moerman (kevinmoerman@hotmail.com)