Added beam twist based on deformation gradient tensor

docs/DEMO_febio_0088_dogbone_sample_01.m

@@ -1,7 +1,7 @@
-%% DEMO_febio_0001_cube_uniaxial
+%% DEMO_febio_0088_dogbone_sample_01.m
 % Below is a demonstration for:
 % 
-% * Building geometry for a cube with hexahedral elements
+% * Building geometry for a dogbone sample 
 % * Defining the boundary conditions 
 % * Coding the febio structure
 % * Running the model

docs/DEMO_febio_0091_beam_twist_F.m

@@ -0,0 +1,375 @@
+%% DEMO_febio_0091_beam_twist_F
+% Below is a demonstration for:
+% 
+% * Building geometry for a beam with hexahedral elements
+% * Defining the boundary conditions 
+% * Coding the febio structure
+% * Running the model
+% * Importing and visualizing the displacement results
+
+%% Keywords
+%
+% * febio_spec version 4.0
+% * febio, FEBio
+% * beam torsion twist loading
+% * prescribed displacement boundary condition
+% * hexahedral elements, hex8, hex20
+% * beam, rectangular
+% * static, solid, multi-step
+% * hyperelastic, Ogden
+% * displacement logfile
+% * stress logfile
+
+%%
+
+clear; close all; clc;
+
+%% Plot settings
+fontSize=20;
+faceAlpha1=0.8;
+markerSize=40;
+markerSize2=20;
+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=[febioFebFileNamePart,'.txt']; %FEBio log file name
+febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
+
+%Specifying dimensions and number of elements
+beamWidth=10; 
+sampleWidth=beamWidth; %Width 
+sampleThickness=4*beamWidth; %Thickness 
+sampleHeight=beamWidth; %Height
+pointSpacings=3*ones(1,3); %Desired point spacing between nodes
+numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
+numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
+numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3
+
+elementType='hex20'; %'hex8'
+
+%Define applied torsion angle
+alphaRotTotal=2*pi;%0.5*pi; %Total twist angle
+numSteps=50; %Number of steps
+
+%Material parameter set
+c1=1e-3; %Shear-modulus-like parameter
+m1=2; %Material parameter setting degree of non-linearity
+k_factor=1e2; %Bulk modulus factor 
+k=c1*k_factor; %Bulk modulus
+
+% FEA control settings
+numTimeSteps=30; %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
+symmetric_stiffness=1; 
+
+runMode='internal';% 'internal' or 'external'
+
+%% Creating model geometry and mesh
+% A box is created with tri-linear hexahedral (hex8) elements using the
+% |hexMeshBox| function. The function offers the boundary faces with
+% seperate labels for the top, bottom, left, right, front, and back sides.
+% As such these can be used to define boundary conditions on the exterior. 
+
+% Create a box with hexahedral elements
+beamDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
+beamElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
+outputStructType=2; %A structure compatible with mesh view
+[meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);
+
+%Access elements, nodes, and faces from the structure
+E=meshStruct.elements; %The elements 
+V=meshStruct.nodes; %The nodes (vertices)
+Fb=meshStruct.facesBoundary; %The boundary faces
+Cb=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
+elementMaterialIndices=ones(size(E,1),1); %Element material indices
+
+if strcmp(elementType,'hex20')
+    [E,V,~,Fb]=hex8_hex20(E,V,{},Fb);
+    meshStruct.elements=E;
+    meshStruct.nodes=V;
+    meshStruct.Fb=Fb;
+    quadType='quad8';
+else
+    quadType='quad4';
+end
+
+%% 
+% Plotting model boundary surfaces and a cut view
+
+hFig=cFigure; 
+
+subplot(1,2,1); hold on; 
+title('Model boundary surfaces and labels','FontSize',fontSize);
+hp=gpatch(Fb,V,Cb,'k',faceAlpha1); 
+hp.Marker='.';
+hp.MarkerSize=markerSize2;
+
+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;
+
+%% Create rigid body end plate for rotation constraint
+Fr=Fb(Cb==3,:); %The faces for the beam end
+
+%Nodes part of the rigid body
+indRigid=unique(Fr(:));
+
+%% Defining the boundary conditions
+% The visualization of the model boundary shows colors for each side of the
+% cube. These labels can be used to define boundary conditions. 
+
+%Define supported node set
+bcSupportList=unique(Fb(Cb==4,:)); %Node set part of selected face
+
+%% 
+% 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,'w','k',0.5);
+
+hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
+plotV(V(indRigid,:),'r.','MarkerSize',markerSize);
+legend(hl,{'BC support'});
+
+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='4.0'; 
+
+%Module section
+febio_spec.Module.ATTR.type='solid'; 
+
+%Control section
+febio_spec.Control.analysis='STATIC';
+febio_spec.Control.time_steps=numTimeSteps;
+febio_spec.Control.step_size=1/numTimeSteps;
+febio_spec.Control.solver.max_refs=max_refs;
+febio_spec.Control.solver.qn_method.max_ups=max_ups;
+febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
+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;
+
+%Material section
+materialName1='Material1';
+febio_spec.Material.material{1}.ATTR.name=materialName1;
+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;
+
+%Mesh section
+% -> Nodes
+febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
+febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
+febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates
+
+% -> Elements
+partName1='Part1';
+febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
+febio_spec.Mesh.Elements{1}.ATTR.type=elementType; %Element type 
+febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
+febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix
+
+% -> NodeSets
+nodeSetName1='bcSupportList';
+febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
+febio_spec.Mesh.NodeSet{1}.VAL=mrow(bcSupportList);
+
+nodeSetName2='bcPrescribedDeformation';
+febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
+febio_spec.Mesh.NodeSet{2}.VAL=mrow(indRigid);
+
+%MeshDomains section
+febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
+febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;
+
+%Boundary condition section 
+% -> Fix boundary conditions
+febio_spec.Boundary.bc{1}.ATTR.name='zero_displacement_xyz';
+febio_spec.Boundary.bc{1}.ATTR.type='zero displacement';
+febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
+febio_spec.Boundary.bc{1}.x_dof=1;
+febio_spec.Boundary.bc{1}.y_dof=1;
+febio_spec.Boundary.bc{1}.z_dof=1;
+
+
+Q = euler2DCM([0,0.25*pi,0]);
+
+% febio_spec.Boundary.bc{2}.ATTR.name='prescribed_deformation';
+febio_spec.Boundary.bc{2}.ATTR.type='prescribed deformation';
+febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
+febio_spec.Boundary.bc{2}.scale.ATTR.lc=1;
+febio_spec.Boundary.bc{2}.scale.VAL=1;
+febio_spec.Boundary.bc{2}.F=Q(:)';
+
+% <bc type="prescribed deformation" node_set="set1">
+%   <scale lc="1">0.1</scale>
+%   <F>2,0,0, 0,2,0, 0,0.25,0</F>
+% <bc>
+
+%LoadData section
+% -> load_controller
+febio_spec.LoadData.load_controller{1}.ATTR.name='LC_1';
+febio_spec.LoadData.load_controller{1}.ATTR.id=1;
+febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
+febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
+%febio_spec.LoadData.load_controller{1}.extend='CONSTANT';
+febio_spec.LoadData.load_controller{1}.points.pt.VAL=[0 0; 1 1];
+
+%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=',';
+
+%% 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
+
+%% 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.runMode=runMode;
+
+[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!
+
+%% Import FEBio results 
+
+if runFlag==1 %i.e. a succesful run
+
+    %% 
+    % Importing nodal displacements from a log file
+    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),0,1);
+
+    %Access data
+    N_disp_mat=dataStruct.data; %Displacement
+    timeVec=dataStruct.time; %Time
+
+    %Create deformed coordinate set
+    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
+
+    %% 
+    % Plotting the simulated results using |anim8| to visualize and animate
+    % deformations 
+
+    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude
+
+    % Create basic view and store graphics handle to initiate animation
+    hf=cFigure; %Open figure  
+    gtitle([febioFebFileNamePart,': Press play to animate']);
+    title('Displacement magnitude [mm]','Interpreter','Latex')
+    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
+    hp.Marker='.';
+    hp.MarkerSize=markerSize2;
+    hp.FaceColor='interp';
+    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object
+
+    axisGeom(gca,fontSize); 
+    colormap(gjet(250)); colorbar;
+    caxis([0 max(DN_magnitude)]);    
+    axis(axisLim(V_DEF)); %Set axis limits statically
+    camlight headlight;        
+
+    % Set up animation features
+    animStruct.Time=timeVec; %The time vector    
+    for qt=1:1:size(N_disp_mat,3) %Loop over time increments        
+        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude
+
+        %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(:,:,qt),DN_magnitude}; %Property values for to set in order to animate
+    end        
+    anim8(hf,animStruct); %Initiate animation feature    
+    drawnow;
+
+end
+
+%% 
+%
+% <<gibbVerySmall.gif>>
+% 
+% _*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) 2006-2023 Kevin Mattheus Moerman and the GIBBON contributors
+% 
+% 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/>.

docs/Gemfile

@@ -1,3 +1,4 @@
 source 'https://rubygems.org'
 gem 'github-pages', group: :jekyll_plugins
 gem "webrick"
+gem "ffi", github: "ffi/ffi"