DEMO_febio_0012_disc_pressure

Below is a demonstration for:

Contents

Keywords

clear; close all; clc;

Plot settings

fontSize=20;
faceAlpha1=0.8;
markerSize=40;
lineWidth=3;

Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting force

%Load
pressureValue=0.5e-3;
loadType='pressure';%'traction';

%Specifying dimensions and number of elements
pointSpacing=1;
inputStruct.cylRadius=30;
inputStruct.numRadial=round((inputStruct.cylRadius*pi)/pointSpacing);
inputStruct.cylHeight=3;
inputStruct.numHeight=round(inputStruct.cylHeight/pointSpacing);
inputStruct.meshType='tri';
inputStruct.closeOpt=1;

% Derive patch data for a cylinder
[Fs,Vs,Cs]=patchcylinder(inputStruct);

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=8; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

Creating model geometry and mesh

The disc is meshed using tetrahedral elements using tetgen

inputStruct.stringOpt='-pq1.2AaY';
inputStruct.Faces=Fs;
inputStruct.Nodes=Vs;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=Cs; %Face boundary markers
inputStruct.regionPoints=getInnerPoint(Fs,Vs); %region points
inputStruct.regionA=tetVolMeanEst(Fs,Vs); %Volume for regular tets
inputStruct.minRegionMarker=2; %Minimum region marker
[meshStruct]=runTetGen(inputStruct); % Mesh model using tetrahedral elements using tetGen
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- TETGEN Tetrahedral meshing --- 04-Jun-2019 12:46:30
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Writing SMESH file --- 04-Jun-2019 12:46:30
----> Adding node field
----> Adding facet field
----> Adding holes specification
----> Adding region specification
--- Done --- 04-Jun-2019 12:46:30
--- Running TetGen to mesh input boundary--- 04-Jun-2019 12:46:30
Opening /mnt/data/MATLAB/GIBBON/data/temp/temp.smesh.
Delaunizing vertices...
Delaunay seconds:  0.035627
Creating surface mesh ...
Surface mesh seconds:  0.004654
Recovering boundaries...
Boundary recovery seconds:  0.00625
Removing exterior tetrahedra ...
Spreading region attributes.
Exterior tets removal seconds:  0.000893
Recovering Delaunayness...
Delaunay recovery seconds:  0.015698
Refining mesh...
Refinement seconds:  0.029619
Optimizing mesh...
Optimization seconds:  0.0019

Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.node.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.ele.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.face.
Writing /mnt/data/MATLAB/GIBBON/data/temp/temp.1.edge.

Output seconds:  0.029795
Total running seconds:  0.124624

Statistics:

  Input points: 1836
  Input facets: 3668
  Input segments: 5502
  Input holes: 0
  Input regions: 1

  Mesh points: 2883
  Mesh tetrahedra: 11786
  Mesh faces: 25406
  Mesh faces on exterior boundary: 3668
  Mesh faces on input facets: 3668
  Mesh edges on input segments: 5502
  Steiner points inside domain: 1047

--- Done --- 04-Jun-2019 12:46:30
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- Importing TetGen files --- 04-Jun-2019 12:46:30
--- Done --- 04-Jun-2019 12:46:30

Access model element and patch data

Fb=meshStruct.facesBoundary;
Cb=meshStruct.boundaryMarker;
V=meshStruct.nodes;
CE=meshStruct.elementMaterialID;
E=meshStruct.elements;

Plotting model boundary surfaces and a cut view

hFig=cFigure;
subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb,V,Cb,'k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);
drawnow;

Defining the boundary conditions

The visualization of the model boundary shows colors for each side of the disc. These labels can be used to define boundary conditions.

%Define supported node sets
logicFace=Cb==1; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==2; %Logic for current face set
F_pressure=fliplr(Fb(logicFace,:)); %The current face set

Visualizing boundary conditions. Markers plotted on the semi-transparent model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','none',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(2)=gpatch(F_pressure,V,'r','k',1);
patchNormPlot(F_pressure,V);
legend(hl,{'BC full support','Pressure surface'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

Defining the FEBio input structure

See also febioStructTemplate and febioStruct2xml and the FEBio user manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='2.5';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis.ATTR.type='static';
febio_spec.Control.title='Disc analysis';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;
febio_spec.Control.max_refs=max_refs;
febio_spec.Control.max_ups=max_ups;

%Material section
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

%Geometry section
% -> Nodes
febio_spec.Geometry.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Geometry.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Geometry.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
febio_spec.Geometry.Elements{1}.ATTR.type='tet4'; %Element type of this set
febio_spec.Geometry.Elements{1}.ATTR.mat=1; %material index for this set
febio_spec.Geometry.Elements{1}.ATTR.name='Disc'; %Name of the element set
febio_spec.Geometry.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Geometry.Elements{1}.elem.VAL=E;

% -> NodeSets
febio_spec.Geometry.NodeSet{1}.ATTR.name='bcSupportList';
febio_spec.Geometry.NodeSet{1}.node.ATTR.id=bcSupportList(:);

% -> Surfaces
febio_spec.Geometry.Surface{1}.ATTR.name='Pressure_surface';
febio_spec.Geometry.Surface{1}.tri3.ATTR.lid=(1:size(F_pressure,1))';
febio_spec.Geometry.Surface{1}.tri3.VAL=F_pressure;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.fix{1}.ATTR.bc='x';
febio_spec.Boundary.fix{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Boundary.fix{2}.ATTR.bc='y';
febio_spec.Boundary.fix{2}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Boundary.fix{3}.ATTR.bc='z';
febio_spec.Boundary.fix{3}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;

%Loads
switch loadType
    case 'pressure'
        febio_spec.Loads.surface_load{1}.ATTR.type='pressure';
        febio_spec.Loads.surface_load{1}.ATTR.surface=febio_spec.Geometry.Surface{1}.ATTR.name;
        febio_spec.Loads.surface_load{1}.pressure.VAL=pressureValue;
        febio_spec.Loads.surface_load{1}.pressure.ATTR.lc=1;
    case 'traction'
        febio_spec.Loads.surface_load{1}.ATTR.type='traction';
        febio_spec.Loads.surface_load{1}.ATTR.surface=febio_spec.Geometry.Surface{1}.ATTR.name;
        febio_spec.Loads.surface_load{1}.scale.VAL=pressureValue;
        febio_spec.Loads.surface_load{1}.scale.ATTR.lc=1;
        febio_spec.Loads.surface_load{1}.traction=[0 0 -1];
end

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

Quick viewing of the FEBio input file structure

The febView function can be used to view the xml structure in a MATLAB figure window.

febView(febio_spec); %Viewing the febio file

Exporting the FEBio input file

Exporting the febio_spec structure to an FEBio input file is done using the febioStruct2xml function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode
% febView(febioFebFileName);

Running the FEBio analysis

To run the analysis defined by the created FEBio input file the runMonitorFEBio function is used. The input for this function is a structure defining job settings e.g. the FEBio input file name. The optional output runFlag informs the user if the analysis was run succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode='external';%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=3; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--- STARTING FEBIO JOB --- 04-Jun-2019 12:46:35
Waiting for log file...
Proceeding to check log file...04-Jun-2019 12:46:37
------- converged at time : 0.1
------- converged at time : 0.181833
------- converged at time : 0.2654
------- converged at time : 0.352254
------- converged at time : 0.441737
------- converged at time : 0.533323
------- converged at time : 0.626592
------- converged at time : 0.721208
------- converged at time : 0.8169
------- converged at time : 0.913454
------- converged at time : 1
--- Done --- 04-Jun-2019 12:46:40

Import FEBio results

if runFlag==1 %i.e. a succesful run
    % Importing nodal displacements from a log file
    [~, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN(:,3).^2,2));
    V_def=V+DN;
    [CF]=vertexToFaceMeasure(Fb,DN_magnitude);

    % Importing element stress from a log file
    [time_mat, E_stress_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress)); %Nodal forces
    time_mat=[0; time_mat(:)]; %Time
    stress_cauchy_sim=[0; mean(squeeze(E_stress_mat(:,end,:)),1)'];

Plotting the simulated results using anim8 to visualize and animate deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(Fb,V_def,CF,'k',1); %Add graphics object to animate
    gpatch(Fb,V,'kw','none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]);
    axis([min(V_def(:,1)) max(V_def(:,1)) min(V_def(:,2)) max(V_def(:,2)) min(V_def(:,3)) max(V_def(:,3))]); %Set axis limits statically
    view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude
        V_def=V+DN; %Current nodal coordinates
        [CF]=vertexToFaceMeasure(Fb,DN_magnitude); %Current color data to use

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,CF}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;
end

GIBBON www.gibboncode.org

Kevin Mattheus Moerman, [email protected]

GIBBON footer text

License: https://github.com/gibbonCode/GIBBON/blob/master/LICENSE

GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for image segmentation, image-based modeling, meshing, and finite element analysis.

Copyright (C) 2019 Kevin Mattheus Moerman

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.