doubleContraction
Below is a demonstration of the features of the doubleContraction function
Contents
clear; close all; clc;
Syntax
C=doubleContraction(A,B,subContract);
Description
This function performs a double contraction of two tensors A and B to produce the output tensor C. The tensors are either both 3x3 second order tensors or A is a 4th order 3x3x3x3 tensor. The optional 3rd input subContract defines the indices to contract over for 4th order tensors e.g. subContract=[1 2] would contract over the first 2 indices.
Examples
Double contraction of a 2nd-order and a 2nd-order tensor
Creating an example tensor. Here a tensor A with all zeros except for the 3,3 component is created. Next a tensor B which, through double contraction, can "pick out" this scalar value is created. Finally the double contraction is performed illustrating this process.
%Create example tensor A A=zeros(3,3); A(3,3)=pi; %Zeros with Azz=A33=pi R=euler2DCM([0.25*pi 0.25*pi 0]); %Rotation matrix A=R*A*R' %Rotated version of A % Create 'structure' tensor example B b=[0 0 1]*R.'; %Rotated Z vector B=b.'*b %Rotated "structure" tensor for vector b
A = 1.5708 -1.1107 1.1107 -1.1107 0.7854 -0.7854 1.1107 -0.7854 0.7854 B = 0.5000 -0.3536 0.3536 -0.3536 0.2500 -0.2500 0.3536 -0.2500 0.2500
Double contraction using doubleContraction
s=doubleContraction(A,B)
s = 3.1416
Double contraction of a 4th-order and a 2nd-order tensor
Example for Hooke's law of linear elasticity
%Creating the 4th order linear elastic stiffness tensor I=eye(3,3); %The 2nd order identity tensor II1=dyadicProduct(I,I,1); %4th order base tensor 1 II3=dyadicProduct(I,I,3); %4th order base tensor 3 %Parameters for Hooke's law mu=1; %The shear modulus lambda=5; %The lambda lame parameter %Compose stiffness tensor C=lambda.*II1+2.*mu.*II3; %Create strain tensor E=rand(3,3); E=E*E'
E = 1.0793 0.8967 1.4848 0.8967 0.8533 1.4394 1.4848 1.4394 2.4337
Double contraction using doubleContraction
subContract=[3 4]; S=doubleContraction(C,E,subContract)
S = 23.9903 1.7933 2.9697 1.7933 23.5381 2.8788 2.9697 2.8788 26.6990
Using symbolic variables
Example for Hooke's law of linear elasticity
%Creating the 4th order linear elastic stiffness tensor I=eye(3,3); %The 2nd order identity tensor II1=dyadicProduct(I,I,1); %4th order base tensor 1 II3=dyadicProduct(I,I,3); %4th order base tensor 3 %Parameters for Hooke's law try syms mu lambda catch mu=1; %The shear modulus lambda=5; %The lambda lame parameter end %Compose stiffness tensor C=lambda.*II1+2.*mu.*II3; %Display Voigt mapped version c=voigtMap(C) %Create strain tensor syms E1_1 E2_2 E3_3 E1_2 E1_3 E2_3; E=[E1_1 E1_2 E1_3;... E1_2 E2_2 E2_3;... E1_3 E2_3 E3_3]
c = [ lambda + 2*mu, lambda, lambda, 0, 0, 0] [ lambda, lambda + 2*mu, lambda, 0, 0, 0] [ lambda, lambda, lambda + 2*mu, 0, 0, 0] [ 0, 0, 0, mu, 0, 0] [ 0, 0, 0, 0, mu, 0] [ 0, 0, 0, 0, 0, mu] E = [ E1_1, E1_2, E1_3] [ E1_2, E2_2, E2_3] [ E1_3, E2_3, E3_3]
Double contraction using doubleContraction
subContract=[3 4]; S=doubleContraction(C,E,subContract)
S = [ E2_2*lambda + E3_3*lambda + E1_1*(lambda + 2*mu), 2*E1_2*mu, 2*E1_3*mu] [ 2*E1_2*mu, E1_1*lambda + E3_3*lambda + E2_2*(lambda + 2*mu), 2*E2_3*mu] [ 2*E1_3*mu, 2*E2_3*mu, E1_1*lambda + E2_2*lambda + E3_3*(lambda + 2*mu)]
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/.