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)