Skip to content

Commit

Permalink
project initialized
Browse files Browse the repository at this point in the history
  • Loading branch information
@MarlonFranke
MarlonFranke committed Sep 29, 2022
1 parent de2b6dc commit 7af299d
Show file tree
Hide file tree
Showing 7,910 changed files with 334,856 additions and 0 deletions.
The diff you're trying to view is too large. We only load the first 3000 changed files.
2 changes: 2 additions & 0 deletions MoofeKIT.prj
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<MATLABProject xmlns="http://www.mathworks.com/MATLABProjectFile" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0"/>
120 changes: 120 additions & 0 deletions commonFunctions/@spectralDecomp/assemMatTang.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
function out=assemMatTang(obj,stress,eMat,factStress)
%% Assembly of material tangent tensor
% member function of spectralDecomp
% always in voigt notation
% Robin Pfefferkorn 23.9.2019

%changelog
%----------------------------------------------------------
%v1.0 (23.9.2019)
%first release
%
%v1.1 (21.10.20)
%added version without perturbation using exact eigen vectors and
%adjusted tangent expressions in case of identical eigenvalues

%% input
if nargin<=3
factStress = -2;
end


%% computations
lam = obj.principalStretch;

%bases
basesVoigt = obj.eigenBasesVoigt;
M11_v = basesVoigt(:,1);
M22_v = basesVoigt(:,2);
M33_v = basesVoigt(:,3);
M12_v = basesVoigt(:,4);
M23_v = basesVoigt(:,5);
M13_v = basesVoigt(:,6);

%% principal stress contributions
if obj.spatialFormulation
principalStr = factStress*stress;
else
principalStr = factStress*(lam.^-2).*stress;
end

%% geometric contributions
switch obj.algorithmNr
%------------------------------------------------------------------
%Perturbation
%------------------------------------------------------------------
case 0
if obj.spatialFormulation
%geometric parts
geom = zeros(3,1);
geom(1) = (stress(2)*lam(1)^2-stress(1)*lam(2)^2)/(lam(2)^2-lam(1)^2);
geom(2) = (stress(3)*lam(2)^2-stress(2)*lam(3)^2)/(lam(3)^2-lam(2)^2);
geom(3) = (stress(3)*lam(1)^2-stress(1)*lam(3)^2)/(lam(3)^2-lam(1)^2);
else
%geometric parts
geom = zeros(3,1);
geom(1) = (stress(2)-stress(1))/(lam(2)^2-lam(1)^2);
geom(2) = (stress(3)-stress(2))/(lam(3)^2-lam(2)^2);
geom(3) = (stress(3)-stress(1))/(lam(3)^2-lam(1)^2);
end
%------------------------------------------------------------------
%exact eigenvectors
%------------------------------------------------------------------
case 1
threshold = obj.perturbThreshold;
if obj.spatialFormulation
%geometric parts
geom = zeros(3,1);
geom(1) = switchGeomPartSpatial(lam,stress,eMat,1,2,threshold);
geom(2) = switchGeomPartSpatial(lam,stress,eMat,2,3,threshold);
geom(3) = switchGeomPartSpatial(lam,stress,eMat,1,3,threshold);
else
%geometric parts
geom = zeros(3,1);
geom(1) = switchGeomPartMaterial(lam,stress,eMat,1,2,threshold);
geom(2) = switchGeomPartMaterial(lam,stress,eMat,2,3,threshold);
geom(3) = switchGeomPartMaterial(lam,stress,eMat,1,3,threshold);
end
otherwise
error('algorithm not implemented')
end

%% assembly
switch obj.algorithmNr
%------------------------------------------------------------------
%perturbation, exact eigenvectors (eigenvectors known)
%------------------------------------------------------------------
case {0,1}
out = ...
eMat(1,1)*(M11_v*M11_v') + eMat(1,2)*(M11_v*M22_v') + eMat(1,3)*(M11_v*M33_v') + ...
eMat(2,1)*(M22_v*M11_v') + eMat(2,2)*(M22_v*M22_v') + eMat(2,3)*(M22_v*M33_v') + ...
eMat(3,1)*(M33_v*M11_v') + eMat(3,2)*(M33_v*M22_v') + eMat(3,3)*(M33_v*M33_v') + ...
principalStr(1)*(M11_v*M11_v')+...
principalStr(2)*(M22_v*M22_v')+...
principalStr(3)*(M33_v*M33_v')+...
geom(1)*(M12_v*M12_v')+...
geom(2)*(M23_v*M23_v')+...
geom(3)*(M13_v*M13_v');
otherwise
error('algorithm not implemented')
end

end
function out = switchGeomPartSpatial(lam,stress,eMat,a,b,threshold)
if abs(lam(a)-lam(b))/lam(a)>threshold
%distinct eigenvalues
out = (stress(b)*lam(a)^2-stress(a)*lam(b)^2)/(lam(b)^2-lam(a)^2);
else
%coincident eigenvalues
out = 1/2*(eMat(b,b)-eMat(a,b))-stress(a);
end
end
function out = switchGeomPartMaterial(lam,stress,eMat,a,b,threshold)
if abs(lam(a)-lam(b))/lam(a)>threshold
%distinct eigenvalues
out = (stress(b)-stress(a))/(lam(b)^2-lam(a)^2);
else
%coincident eigenvalues
out = 1/2*(eMat(b,b)-eMat(a,b));
end
end
28 changes: 28 additions & 0 deletions commonFunctions/@spectralDecomp/assemTens.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
function out=assemTens(obj,stress)
%% Assembly of second order tensor
% member function of spectralDecomp
% Robin Pfefferkorn 23.9.2019

%changelog
%----------------------------------------------------------
%v1.0 (23.9.2019)
%first release
%
%v1.1 (21.10.20)
%added version without perturbation using exact eigen vectors and
%adjusted tangent expressions in case of identical eigenvalues


%% computations
switch obj.algorithmNr
%------------------------------------------------------------------
%perturbation, exact eigenvectors
%------------------------------------------------------------------
case {0,1}
eigenBases = obj.eigenBases;
out = stress(1)*eigenBases{1} + stress(2)*eigenBases{2} + stress(3)*eigenBases{3};
otherwise
error('algorithm not implemented')
end

end
109 changes: 109 additions & 0 deletions commonFunctions/@spectralDecomp/decomposition.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
function obj=decomposition(obj,b,mapVoigt)
%% Spectral Decomposition Algorithms
% member function of spectralDecomp
% Robin Pfefferkorn 23.9.2019
%
% INPUT:
% b...right or left cauchy green tensor (depending on spatial or
% material formulation)
% mapVoigt...vector containing voigt map information
%
%changelog
%----------------------------------------------------------
%v1.0 (23.9.2019)
%first release
%
%v1.1 (21.10.20)
%added version without perturbation using exact eigen vectors and
%adjusted tangent expressions in case of identical eigenvalues

%% check input
if size(b,1)~=3 || size(b,2)~=3
error('only 3x3 matrices supported')
end

temp=nargin;
if temp==3
voigt=true;
if isrow(mapVoigt)
mapVoigt = mapVoigt';
elseif ~iscolumn(mapVoigt)
error('mapVoigt must be a column vector')
end
elseif temp==2
voigt=false;
end

%% compute eigenvalues and vectors
switch obj.algorithmNr
%------------------------------------------------------------------
%Perturbation
%------------------------------------------------------------------
case 0
%matlab decomp
[nEig,lam2] = eig(b);
lam = sqrt(diag(lam2));

%perturbations
tol = obj.perturbThreshold*max(abs(lam));
pertubed = tol;
if abs(lam(1)-lam(2)) < tol
lam(1) = (1-pertubed)*lam(1);
lam(2) = (1+pertubed)*lam(2);
lam(3) = 1/((1-pertubed)*(1+pertubed))*lam(3);
elseif abs(lam(1)-lam(3)) < tol
lam(1) = (1-pertubed)*lam(1);
lam(2) = 1/((1-pertubed)*(1+pertubed))*lam(2);
lam(3) = (1+pertubed)*lam(3);
elseif abs(lam(2)-lam(3)) < tol
lam(1) = 1/((1-pertubed)*(1+pertubed))*lam(1);
lam(2) = (1-pertubed)*lam(2);
lam(3) = (1+pertubed)*lam(3);
end
%------------------------------------------------------------------
%exact eigenvectors
%------------------------------------------------------------------
case 1
%matlab decomp
[nEig,lam2] = eig(b);
lam = sqrt(diag(lam2));
otherwise
error('algorithm not implemented')
end

%% compute bases and save results
switch obj.algorithmNr
%------------------------------------------------------------------
%perturbation, exact eigenvectors (eigenvectors known)
%------------------------------------------------------------------
case {0,1}
%save results
obj.principalStretch = lam;
obj.eigenVectors = nEig;

%bases
bases = cell(6,1);
bases{1} = nEig(:,1)*nEig(:,1).';
bases{2} = nEig(:,2)*nEig(:,2).';
bases{3} = nEig(:,3)*nEig(:,3).';
bases{4} = nEig(:,1)*nEig(:,2).'+nEig(:,2)*nEig(:,1).';
bases{5} = nEig(:,2)*nEig(:,3).'+nEig(:,3)*nEig(:,2).';
bases{6} = nEig(:,1)*nEig(:,3).'+nEig(:,3)*nEig(:,1).';
obj.eigenBases = bases;

%bases - voigt
if voigt
basesVoigt = zeros(numel(mapVoigt),6);
basesVoigt(:,1) = bases{1}(mapVoigt);
basesVoigt(:,2) = bases{2}(mapVoigt);
basesVoigt(:,3) = bases{3}(mapVoigt);
basesVoigt(:,4) = bases{4}(mapVoigt);
basesVoigt(:,5) = bases{5}(mapVoigt);
basesVoigt(:,6) = bases{6}(mapVoigt);
obj.eigenBasesVoigt = basesVoigt;
end
otherwise
error('algorithm not implemented')
end

end
36 changes: 36 additions & 0 deletions commonFunctions/@spectralDecomp/scaleEigenvectors.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
function obj=scaleEigenvectors(obj,factors)
%% Scaling of eigenvectors with given factors
% member function of spectralDecomp
% Robin Pfefferkorn 23.9.2019

%changelog
%----------------------------------------------------------
%v1.0 (23.7.2020)
%first release
%
%v1.1.....


%% checks
if numel(factors)~=3 || ~isfloat(factors)
error('input must have 3 float entries')
end

% if any(abs(factors-1)>0.01)
% keyboard
% end

%% computations
obj.eigenVectors = obj.eigenVectors*diag(factors);

% bases
factBases = [factors(1)^2, factors(2)^2, factors(3)^2, ...
factors(1)*factors(2), factors(2)*factors(3), factors(1)*factors(3)];

bases = obj.eigenBases;
for i=1:numel(bases)
bases{i} = bases{i}*factBases(i);
end
obj.eigenBases = bases;
obj.eigenBasesVoigt = obj.eigenBasesVoigt*diag(factBases);
end
Loading

0 comments on commit 7af299d

Please sign in to comment.