%This MATLAB-script is used to calculate the end-diastole volume,
%end-systole volume and ejection fraction for the left ventricle model.
%The scrpit uses an Abaqus inputfile which must contain a node set with all
%nodes on endocaridal surface.
%Further at txt file (created with the python script) containing the node
%coordinates at the end-systole state.
clear all close all clc
%Extract relevant line numbers from input file. Change keywords to coincide
counter = counter + 1;
if ischar(tline)
U = strfind(tline, ’*Node’);
if isfinite(U) == 1;
nodeLine = counter;
end
U = strfind(tline, ’*Element’);
if isfinite(U) == 1;
nnodes=elsLine-nodeLine-1;
%Extract node coordinates at end-diastole state.
orgNode = csvread(’TEM_8layers_fnsn150_f45.inp’,nodeLine, 1, [ nodeLine,1,nodeLine+nnodes-1,3]);
%Rotate nodes 90 degrees about x-axis as done with the instance in the
%Abaqus model
tempY = orgNode(:,2);
tempZ = orgNode(:,3);
orgNode(:,2) = cos(pi()/2)*tempY(:)-sin(pi()/2)*tempZ(:);
orgNode(:,3) = sin(pi()/2)*tempY(:)+cos(pi()/2)*tempZ(:);
%Extract node numbers from the endoNode set
endoNodeSet = csvread(’TEM_8layers_fnsn150_f45.inp’,endoNodeline, 0, [endoNodeline,0,endoNodestop-2,15]);
endoNodeSet(endoNodeSet==0) = [];
endoNodeSet = sort(endoNodeSet(:));
%Extract node coordinates at end-systle state.
deformedNode =textread(’deformed_shape_TEM_8layers_fnsn150_f45.txt
’,’’,’headerlines’, 3);
%Creating variables containing node coordinates only from the nodes on the
%Caclutate the end-diastole and end-systole volume using the convhull
%Write out the volumes (in ml) and ejection fraction V*0.001
V2*0.001
ejectionFraction=(V-V2)/V
%The remainder of the script is plots used for visual confirmation of the
%validity of the result figure()
hold on
scatter3(orgEndoCoor(:,1), orgEndoCoor(:,2), orgEndoCoor(:,3), ’.’
)
scatter3(deformedEndoCoor(:,1), deformedEndoCoor(:,2), deformedEndoCoor(:,3), ’r.’)
xlabel(’X’) ylabel(’Y’) zlabel(’Z’) hold off figure() hold on
scatter3(deformedEndoCoor(:,1), deformedEndoCoor(:,2), deformedEndoCoor(:,3), ’.’)
xlabel(’X’) ylabel(’Y’) zlabel(’Z’) hold off figure() hold on
trisurf(TriIdx, orgEndoCoor(:,1),orgEndoCoor(:,2),orgEndoCoor(:,3) )
xlabel(’X’) ylabel(’Y’) zlabel(’Z’) hold off figure()
trisurf(TriIdx2, deformedEndoCoor(:,1),deformedEndoCoor(:,2), deformedEndoCoor(:,3))
xlabel(’X’) ylabel(’Y’) zlabel(’Z’)
ening
#This python script extracts nodal displacement and calculates the longitudinal shortening, radial shortening and wall
thickening
#Variables which needs to be changed are: odb file-path, instance name and node numbers.
from odbAccess import * from abaqusConstants import *
# create odb object from odb file
outputDatabase = openOdb(path=’C:\Users\Gaute\Documents\
MASTERTHESIS\Abaqus\MainModel\TEM_4layers_f100_f60_14.odb’)
# get access to the nodal displacement data
frame = outputDatabase.steps[ ’static’ ].frames[-1]
dispField = frame.fieldOutputs[’U’]
# get access to the part instance
my_part_instance = outputDatabase.rootAssembly.instances[’TEM-4 LAYERS3-1’]
#Variables with node numbers on the form [endocardium,epicardium]
equNodes = [36,27]
ApexNodes = [11, 21]
#Extracting the coordinates before and after deformation
#Note that the note label is not the same as the node number.
numNodesTotal = len( my_part_instance.nodes ) for i in range( numNodesTotal ):
nodeLabel = my_part_instance.nodes[i].label if nodeLabel==equNodes[0]:
node1int = i
node1Pos = my_part_instance.nodes[i].coordinates defNode1Pos = node1Pos + dispField.values[i].data if nodeLabel==equNodes[1]:
node2int = i
node2Pos = my_part_instance.nodes[i].coordinates defNode2Pos = node2Pos + dispField.values[i].data if nodeLabel==ApexNodes[0]:
apexint = i
apexEndoPos = my_part_instance.nodes[i].
coordinates
defEndoApexPos = apexEndoPos + dispField.values[i ].data
if nodeLabel==ApexNodes[1]:
apexEpiint = i
apexEpiPos = my_part_instance.nodes[i].coordinates defEpiApexPos = apexEpiPos + dispField.values[i].
data
#Calculate the wanted parameters
initalWallTh = sqrt((node2Pos[0]-node1Pos[0])**2 + (node2Pos[1]-node1Pos[1])**2 + (node2Pos[2]-node1Pos[2])**2)
defWallTh = sqrt((defNode2Pos[0]-defNode1Pos[0])**2 + (defNode2Pos [1]-defNode1Pos[1])**2)
thFrac = ( defWallTh - initalWallTh) / initalWallTh
initalApexWallTh = sqrt((apexEndoPos[0]-apexEpiPos[0])**2 + ( apexEndoPos[1]-apexEpiPos[1])**2 + (apexEndoPos[2]-apexEpiPos [2])**2)
defApexWallTh = sqrt((defEndoApexPos[0]-defEpiApexPos[0])**2 + ( defEndoApexPos[1]-defEpiApexPos[1])**2 + (defEndoApexPos[2]-defEpiApexPos[2])**2)
thApexFrac = ( defApexWallTh - initalApexWallTh) / initalApexWallTh
initalLong = 20 - apexEndoPos[2]
defLong = 20 - defEndoApexPos[2]
relLongShort = (initalLong - defLong) / initalLong initalRad = node1Pos[0]
defRad = defNode1Pos[0]
radShort = (initalRad-defRad)/initalRad
#Printing the parameters to the command window print "equ wallth:"
F.1 Python Script
#This Python script extracts displacement history of a set of chosen nodes as part of calculating the rotation of the left ventricle model.
#The script writes the node coordinates for each frame of the analysis to a file. This must be done for the 4 nodes on the endocardium and epicardium of both the base and apex.
#The four files created using this script is processed further using a MATLAB script.
from odbAccess import * from abaqusConstants import * from math import *
# create odb object from odb file
outputDatabase = openOdb(path=’C:\Users\Gaute\Documents\
MASTERTHESIS\Abaqus\MainModel\TEM_4layers_fnsn150_f45_14.odb’)
# get access to the nodal displacement data
allFrames = outputDatabase.steps[ ’static’ ].frames lastFrame = outputDatabase.steps[ ’static’ ].frames[-1]
numFrames = len (allFrames)
dispField = lastFrame.fieldOutputs[’U’]
# get access to the part instance
my_part_instance = outputDatabase.rootAssembly.instances[’TEM-4 LAYER3-1’]
#Declare variables with the node numbers from the different regions on the form [center, left, up, right down ] when looking from the apex to the base.
#The center basal node number is set to 0, and given the global coordinate (0,0,20) which is the center of the base in the model.
apexEpiNodes = [21, 419, 376, 496, 615 ] apexEndoNodes = [11, 553, 144, 187, 345 ] basalEndoNodes = [0, 40, 51, 50, 41 ] basalEpiNodes = [0, 49, 47, 44, 48 ]
#Choose node data to write to file activeNodes = apexEpiNodes8layers
#Write data to files (change file name to coincide with nodes)
outFile = open("rotation_TEM_8layers_fnsn150_f45_apexEpi.txt" , ’w
’ )
# write points
numNodesTotal = len( my_part_instance.nodes ) for i in range( numNodesTotal ):
nodeLabel = my_part_instance.nodes[i].label for k in range (numFrames):
if activeNodes[0]==0:
defNodeCenterPos = [0,0,20]
else:
nodeCenterDisp = outputDatabase.steps[ ’static’ ].
frames[k].fieldOutputs[’U’].values[
nodeCenterint].data
defNodeCenterPos = nodeCenterPos + nodeCenterDisp nodeLeftDisp = outputDatabase.steps[ ’static’ ].frames[k].
fieldOutputs[’U’].values[nodeLeftint].data defNodeLeftPos = nodeLeftPos + nodeLeftDisp
nodeUpDisp = outputDatabase.steps[ ’static’ ].frames[k].
fieldOutputs[’U’].values[nodeUpint].data defNodeUpPos = nodeUpPos + nodeUpDisp
nodeRightDisp = outputDatabase.steps[ ’static’ ].frames[k ].fieldOutputs[’U’].values[nodeRightint].data
defNodeRightPos = nodeRightPos + nodeRightDisp
nodeDownDisp = outputDatabase.steps[ ’static’ ].frames[k].
fieldOutputs[’U’].values[nodeDownint].data
for j in range( 3 ):
outFile.write( str( defNodeCenterPos[j] ) + ’ ’ ) for j in range( 3 ):
outFile.write( str( defNodeLeftPos[j] ) + ’ ’ ) for j in range( 3 ):
outFile.write( str( defNodeUpPos[j] ) + ’ ’ ) for j in range( 3 ):
outFile.write( str( defNodeRightPos[j] ) + ’ ’ ) for j in range( 3 ):
if j<2:
outFile.write( str( defNodeDownPos[j] ) +
’ ’ ) else:
outFile.write( str( defNodeDownPos[j] ) ) outFile.write("\n")
outputDatabase.close()
F.2 MATLAB Script
%This MATLAB script calculates the rotation of endocardium and epicardium
%for both the base and apex for all frames through the Abaqus analysis. It
%takes as input 4 txt files with nodal coordinate data, created with the
%Python script.
%NOTE: the function atan2 change sign at 180 degrees which is relevant for
%some of the nodes. For simplicity, when this is relevant, the node in
%question is not included in the calculation of the average rotation. As
%the rotation of the different nodes are almost equal, this has little or
%no effect on the results.
clear all close all clc
%Reading node data from files nodeCoorApexEpi =textread(’
rotation_TEM_4layers_fnsn150_f45_apexEpi.txt’,’’,’headerlines’
, 0);
angleV=zeros(4); %Angles of the vectors the 4 undeformed nodes creates in the global csys
defv=zeros(4,2); %X and Y components of the deformed vectors angleDefV=zeros(4,length(nodeCoorApexEpi)); %Angles of the vectors
the 4 deformed nodes creates in the global csys, for every frame
rot=zeros(4,length(nodeCoorApexEpi)); %Rotation of the 4 nodes for every frame
%Caclulate rotation of apex at epicardium for i=1:4
angleDefV = angleDefV - 360*sign(angleDefV);
end
rot(i,j) = (angleDefV(i,j)-angleV(i));
end end
%Average rotation of apex at epicaridum
avgRotApexEpi = -(rot(1,:)+rot(2,:)+rot(3,:)+rot(4,:))/4;
%Caclulate rotation of apex at endocardium
(1,1+3*i))*180/pi();
angleDefV = angleDefV - 360*sign(angleDefV);
end
rot(i,j) = (angleDefV(i,j)-angleV(i));
end end
%Average rotation of apex at endocaridum
avgRotApexEndo = -(rot(1,:)+rot(2,:)+rot(4,:))/4;
%Caclulate rotation of base at epicardium for i=1:4
angleDefV = angleDefV - 360*sign(angleDefV);
end
rot(i,j) = (angleDefV(i,j)-angleV(i));
end end
%Average rotation of base at epicaridum
avgRotBasalEpi = -(rot(1,:)+rot(2,:)+rot(3,:)+rot(4,:))/4;
%Caclulate rotation of base at endocardium
angleDefV = angleDefV - 360*sign(angleDefV);
end
rot(i,j) = (angleDefV(i,j)-angleV(i));
end end
%Average rotation of base at endocaridum
avgRotBasalEndo = -(rot(1,:)+rot(2,:)+rot(4,:))/3;
%Write out rotation at end-systole avgRotApexEpi(length(nodeCoorApexEpi)) avgRotBasalEpi(length(nodeCoorBasalEpi)) avgRotApexEndo(length(nodeCoorApexEpi)) avgRotBasalEndo(length(nodeCoorBasalEpi))
%Plot rotatation of all frames for visual confirmation of calculation
legend(’Apex Epi’, ’Apex Endo’, ’Basal Epi’, ’Basal Endo ’) hold off