% Clean slate 
clear
close all
clc
tic

%Threshold the Cube
thickness = 1;
Prob_of_Bridge = 0;
Prob_of_Block = 0.545;
Prob_of_Space = 1-Prob_of_Block;


%Create a block of the walls

Slices = [1 1 0 0 1 1 1 0 1 1 0 0 1 1 0 0 1 1 1 0 1 1];


Cube_size = length(Slices);
Cube_wall = repmat(Slices, Cube_size, 1, Cube_size);
Cube_bridge = repmat(Slices, Cube_size, 1, Cube_size);

%Calculate how many cubes to add
Num_of_spaces = (Cube_size^3)-sum(sum(sum(Cube_wall)));
Num_of_Blocks = round(Prob_of_Block*Num_of_spaces);

Dens_Cube = sum(sum(sum(Cube_wall)))/numel(Cube_wall);
[x_size, y_size, z_size] = size(Cube_wall);
%Loop until density is <= as specified

while Dens_Cube > Prob_of_Block
    %Pick a random element
    %Elem = [randi(x_size),randi(y_size),randi(z_size)];
    %Check that element is solid
    Cube_wall(randi(x_size),randi(y_size),randi(z_size)) = 0;
    %Re-caclulate the density
    Dens_Cube = sum(sum(sum(Cube_wall)))/numel(Cube_wall);
end

%% Add pore channel defects

%Calculate how many cubes to add
Num_of_spaces = (Cube_size^3)-sum(sum(sum(Cube_bridge)));
Num_of_Bridge = round(Prob_of_Bridge*Num_of_spaces);

Block_Count = 0;
while Block_Count < Num_of_Bridge
    Vec = [randi(Cube_size),randi(Cube_size),randi(Cube_size)];
    if Cube_bridge(Vec(1),Vec(2),Vec(3)) == 1
        %You are in a layer, do nothing
    else
        %You are in a gap, evaluate the location
        Space_check = 0;
        if Vec(1) ~= 1
            Space_check = Space_check + Cube_bridge((Vec(1)-1),Vec(2),Vec(3));
        end
        if Vec(1) ~= length(Slices)
            Space_check = Space_check + Cube_bridge((Vec(1)+1),Vec(2),Vec(3));
        end
        if Vec(2) ~= 1
            Space_check = Space_check + Cube_bridge(Vec(1),(Vec(2)-1),Vec(3));
        end
        if Vec(2) ~= length(Slices)
            Space_check = Space_check + Cube_bridge(Vec(1),(Vec(2)+1),Vec(3));
        end        
        if Vec(3) ~= 1
            Space_check = Space_check + Cube_bridge(Vec(1),Vec(2),(Vec(3)-1));
        end
        if Vec(3) ~= length(Slices)
            Space_check = Space_check + Cube_bridge(Vec(1),Vec(2),(Vec(3)+1));
        end
        if Space_check > 0
            Cube_bridge(Vec(1),Vec(2),Vec(3)) = 2;
            Block_Count = Block_Count + 1;
        end
        clear Space_check
    end
end
%%

Cube = Cube_bridge+Cube_wall;
Cube(Cube==1) = 0;
Cube(Cube==2) = 1;

Volume = sum(sum(sum(Cube)))/(Cube_size^3);


%%%%%%%% [BEN's CALCULATE SURFACE AREA CODE] %%%%%%%%%%%%%

Surface_area = 0;
for i = 1:Cube_size
    for j = 1:Cube_size
        for k = 1:Cube_size-1
            if Cube(i,j,k) ~= Cube(i,j,k+1)
                Surface_area = Surface_area + 1;
            end
        end
    end
end
for i = 1:Cube_size
    for j = 1:Cube_size
        for k = 1:Cube_size-1
            if Cube(i,k,j) ~= Cube(i,k+1,j)
                Surface_area = Surface_area + 1;
            end
        end
    end
end
for i = 1:Cube_size
    for j = 1:Cube_size
        for k = 1:Cube_size-1
            if Cube(k,i,j) ~= Cube(k+1,i,j)
                Surface_area = Surface_area + 1;
            end
        end
    end
end

% Varying elements
n = 1;
for i = 1:(Cube_size^3)
    if Cube(i) == 1
        Cube(i) = n;
        n = n+1;
    end
end

maxElem = max(max(max(Cube)));
Cube = Cube+1;


left_string1 = 'C:\Users\vobui\Desktop\Perfect wall\assign';
mid_string = '_w_';
mid_string2 = '_b_';
right_string1 = '.txt';


fid = fopen([left_string1 num2str(thickness) mid_string num2str(Prob_of_Block*1000) mid_string2 num2str(Prob_of_Bridge*10) right_string1],'wt');
    for elem = 1:(Cube_size^3)
        str_1 = 'MPCHG,%1.0f,%1.0f\n';
        fprintf(fid,str_1,Cube(elem),elem);
    end
fclose(fid);
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
left_string2 = 'C:\Users\vobui\Desktop\Perfect wall\materials';
fid = fopen([left_string2 num2str(thickness) mid_string num2str(Prob_of_Block*1000) mid_string2 num2str(Prob_of_Bridge*10) right_string1],'wt');

air_1 = 'MAT,1';
air_2 = 'MP, DENS,,1.3';
air_3 = 'MP, PERX,,1';
air_4 = 'MP, EX,,1E-12';
air_5 = 'MP, PRXY,,1E-12';
fprintf(fid, air_1);
fprintf(fid, '\n');
fprintf(fid, air_2);
fprintf(fid, '\n');
fprintf(fid, air_3);
fprintf(fid, '\n');
fprintf(fid, air_4);
fprintf(fid, '\n');
fprintf(fid, air_5);
fprintf(fid, '\n');


for n = 2:max(max(max(Cube))) % Evaluating material properties
    %% Material properties:
    % Piezoelectric strain d:
    d33 = 443e-12;
    d31 = -208e-12;
    d15 = 669e-12;
    
    % Dielectric permittivity
    perm_x = 906;
    perm_y = 906;
    perm_z = 823;
    
    % Mechanical stiffness
    c11 = 13.4e10;
    c12 = 8.89e10;
    c13 = 9.09e10;
    c33 = 12.1e10;
    c44 = 2.24e10;
    c66 = 2.05e10;
    
    % Assembling matrices & tensors
    d(:,:,1) = [0 0 d15; 0 0 0; d15 0 0];
    d(:,:,2) = [0 0 0; 0 0 d15; 0 d15 0];
    d(:,:,3) = [d31 0 0; 0 d31 0; 0 0 d33];
    c = [c11 c12 c13 0 0 0; c12 c11 c13 0 0 0; c13 c13 c33 0 0 0; 0 0 0 c44 0 0; 0 0 0 0 c44 0; 0 0 0 0 0 c66];


    % Initial random polarisation state
    alpha = 360*rand();
    beta = 360*rand();
    gamma = 0;

    c1 = cosd(alpha); s1 = sind(alpha);
    c2 = cosd(beta); s2 = sind(beta);
    c3 = cosd(gamma); s3 = sind(gamma);
    r = [c1*c2*c3-s1*s3 -c3*s1-c1*c2*s3 c1*s2; c1*s3+c2*c3*s1 c1*c3-c2*s1*s3 s1*s2; -c3*s2 s2*s3 c2]; %Z3 then Y then Z1
    
    %% Matrices transformation
    % Piezo-strain transformation
    d_rot = transform(d, r);
    d_rot_full = [d_rot(1,:,1) d_rot(2,:,1) d_rot(3,:,1); d_rot(1,:,2) d_rot(2,:,2) d_rot(3,:,2); d_rot(1,:,3) d_rot(2,:,3) d_rot(3,:,3)];
    d_rot_reduced = [d_rot_full(:,1) d_rot_full(:,5) d_rot_full(:,9) d_rot_full(:,2) d_rot_full(:,6) d_rot_full(:,3)];
    
    % T_sigma calculation
    Ts1 = [r(1,1)^2 r(2,1)^2 r(3,1)^2 2*r(1,1)*r(2,1) 2*r(2,1)*r(3,1) 2*r(3,1)*r(1,1)];
    Ts2 = [r(1,2)^2 r(2,2)^2 r(3,2)^2 2*r(1,2)*r(2,2) 2*r(2,2)*r(3,2) 2*r(3,2)*r(1,2)];
    Ts3 = [r(1,3)^2 r(2,3)^2 r(3,3)^2 2*r(1,3)*r(2,3) 2*r(2,3)*r(3,3) 2*r(3,3)*r(1,3)];
    Ts4 = [r(1,1)*r(1,2) r(2,1)*r(2,2) r(3,1)*r(3,2) r(1,1)*r(2,2)+r(1,2)*r(2,1) r(2,1)*r(3,2)+r(3,1)*r(2,2) r(3,1)*r(1,2)+r(3,2)*r(1,1)];
    Ts5 = [r(1,2)*r(1,3) r(2,2)*r(2,3) r(3,3)*r(3,2) r(2,3)*r(1,2)+r(1,3)*r(2,2) r(2,2)*r(3,3)+r(3,2)*r(2,3) r(3,2)*r(1,3)+r(1,2)*r(3,3)];
    Ts6 = [r(1,1)*r(1,3) r(2,3)*r(2,1) r(3,3)*r(3,1) r(1,3)*r(2,1)+r(1,1)*r(2,3) r(2,3)*r(3,1)+r(2,1)*r(3,3) r(3,3)*r(1,1)+r(3,1)*r(1,3)];
    Ts = [Ts1;Ts2;Ts3;Ts4;Ts5;Ts6];
    
    % T_epsilon calculation
    Te1 = [r(1,1)^2 r(2,1)^2 r(3,1)^2 r(1,1)*r(2,1) r(2,1)*r(3,1) r(3,1)*r(1,1)];
    Te2 = [r(1,2)^2 r(2,2)^2 r(3,2)^2 r(1,2)*r(2,2) r(2,2)*r(3,2) r(3,2)*r(1,2)];
    Te3 = [r(1,3)^2 r(2,3)^2 r(3,3)^2 r(1,3)*r(2,3) r(2,3)*r(3,3) r(3,3)*r(1,3)];
    Te4 = [2*r(1,1)*r(1,2) 2*r(2,1)*r(2,2) 2*r(3,1)*r(3,2) r(1,1)*r(2,2)+r(1,2)*r(2,1) r(2,1)*r(3,2)+r(3,1)*r(2,2) r(3,1)*r(1,2)+r(3,2)*r(1,1)];
    Te5 = [2*r(1,2)*r(1,3) 2*r(2,2)*r(2,3) 2*r(3,3)*r(3,2) r(2,3)*r(1,2)+r(1,3)*r(2,2) r(2,2)*r(3,3)+r(3,2)*r(2,3) r(3,2)*r(1,3)+r(1,2)*r(3,3)];
    Te6 = [2*r(1,1)*r(1,3) 2*r(2,3)*r(2,1) 2*r(3,3)*r(3,1) r(1,3)*r(2,1)+r(1,1)*r(2,3) r(2,3)*r(3,1)+r(2,1)*r(3,3) r(3,3)*r(1,1)+r(3,1)*r(1,3)];
    Te = [Te1;Te2;Te3;Te4;Te5;Te6];
    
    % Evaluation of rotated stiffness and rotated stress-charge [e] matrix
    c_rot = Ts^-1*c*Te;
    e_rot_ANSYS = transpose(d_rot_reduced*c_rot);
    d_rot_ANSYS = transpose(d_rot_reduced);
    d_rot_ANSYS = [d_rot_ANSYS(1,:);d_rot_ANSYS(2,:);d_rot_ANSYS(3,:);d_rot_ANSYS(5,:);d_rot_ANSYS(4,:);d_rot_ANSYS(6,:)];
    perm = [perm_x; perm_y; perm_z];
    perm_rot = abs(r*perm);
    
    %% Export to APDL commands
    % Permittivity
    p_join_4 = ['MP,PERX,',num2str(n) ',', num2str(perm_rot(1))];
    p_join_5 = ['MP,PERY,',num2str(n) ',', num2str(perm_rot(2))];
    p_join_6 = ['MP,PERZ,',num2str(n) ',', num2str(perm_rot(3))];
    
    element_id_c_poled = ['TB,ANEL,',num2str(n),',1,,0'];
    c_join1 = ['TBDATA,1,',num2str(c_rot(1,1)) ',',num2str(c_rot(1,2)) ',',num2str(c_rot(1,3))];
    c_join2 = ['TBDATA,4,',num2str(c_rot(1,4)) ',',num2str(c_rot(1,5)) ',',num2str(c_rot(1,6))];
    c_join3 = ['TBDATA,7,',num2str(c_rot(2,2)) ',',num2str(c_rot(2,3)) ',',num2str(c_rot(2,4))];
    c_join4 = ['TBDATA,10,',num2str(c_rot(2,5)) ',',num2str(c_rot(2,6)) ',',num2str(c_rot(3,3))];
    c_join5 = ['TBDATA,13,',num2str(c_rot(3,4)) ',',num2str(c_rot(3,5)) ',',num2str(c_rot(3,6))];
    c_join6 = ['TBDATA,16,',num2str(c_rot(4,4)) ',',num2str(c_rot(4,5)) ',',num2str(c_rot(4,6))];
    c_join7 = ['TBDATA,19,',num2str(c_rot(5,5)) ',',num2str(c_rot(5,6)) ',',num2str(c_rot(6,6))];
    
    element_id_e_poled = ['TB,PIEZ,',num2str(n)];
    e_join1 = ['TBDATA,1,',num2str(e_rot_ANSYS(1,1)) ',',num2str(e_rot_ANSYS(1,2)) ',',num2str(e_rot_ANSYS(1,3))];
    e_join2 = ['TBDATA,4,',num2str(e_rot_ANSYS(2,1)) ',',num2str(e_rot_ANSYS(2,2)) ',',num2str(e_rot_ANSYS(2,3))];
    e_join3 = ['TBDATA,7,',num2str(e_rot_ANSYS(3,1)) ',',num2str(e_rot_ANSYS(3,2)) ',',num2str(e_rot_ANSYS(3,3))];
    e_join4 = ['TBDATA,10,',num2str(e_rot_ANSYS(4,1)) ',',num2str(e_rot_ANSYS(4,2)) ',',num2str(e_rot_ANSYS(4,3))];
    e_join5 = ['TBDATA,13,',num2str(e_rot_ANSYS(5,1)) ',',num2str(e_rot_ANSYS(5,2)) ',',num2str(e_rot_ANSYS(5,3))];
    e_join6 = ['TBDATA,16,',num2str(e_rot_ANSYS(6,1)) ',',num2str(e_rot_ANSYS(6,2)) ',',num2str(e_rot_ANSYS(6,3))];
    
    fprintf(fid, '%s', p_join_4);
    fprintf(fid, '\n');
    fprintf(fid, '%s', p_join_5);
    fprintf(fid, '\n');
    fprintf(fid, '%s', p_join_6);
    fprintf(fid, '\n');
    fprintf(fid, '%s', element_id_c_poled);
    fprintf(fid, '\n');
    fprintf(fid, c_join1);
    fprintf(fid, '\n');
    fprintf(fid, c_join2);
    fprintf(fid, '\n');
    fprintf(fid, c_join3);
    fprintf(fid, '\n');
    fprintf(fid, c_join4);
    fprintf(fid, '\n');
    fprintf(fid, c_join5);
    fprintf(fid, '\n');
    fprintf(fid, c_join6);
    fprintf(fid, '\n');
    fprintf(fid, c_join7);
    fprintf(fid, '\n');
    fprintf(fid, '%s', element_id_e_poled);
    fprintf(fid, '\n');
    fprintf(fid, e_join1);
    fprintf(fid, '\n');
    fprintf(fid, e_join2);
    fprintf(fid, '\n');
    fprintf(fid, e_join3);
    fprintf(fid, '\n');
    fprintf(fid, e_join4);
    fprintf(fid, '\n');
    fprintf(fid, e_join5);
    fprintf(fid, '\n');
    fprintf(fid, e_join6);
    fprintf(fid, '\n');
end





for n = 2:max(max(max(Cube))) % Evaluating poled state 
    %% Material properties:

    % Piezoelectric strain d:
    d33 = 443e-12;
    d31 = -208e-12;
    d15 = 669e-12;
    
    % Dielectric permittivity
    perm_x = 906;
    perm_y = 906;
    perm_z = 823;
    
    % Mechanical stiffness
    c11 = 13.4e10;
    c12 = 8.89e10;
    c13 = 9.09e10;
    c33 = 12.1e10;
    c44 = 2.24e10;
    c66 = 2.05e10;
    
    % Assembling matrices & tensors
    d(:,:,1) = [0 0 d15; 0 0 0; d15 0 0];
    d(:,:,2) = [0 0 0; 0 0 d15; 0 d15 0];
    d(:,:,3) = [d31 0 0; 0 d31 0; 0 0 d33];
    c = [c11 c12 c13 0 0 0; c12 c11 c13 0 0 0; c13 c13 c33 0 0 0; 0 0 0 c44 0 0; 0 0 0 0 c44 0; 0 0 0 0 0 c66];


    % Initial random polarisation state
    alpha = 0;
    beta = 0;
    gamma = 0;

    c1 = cosd(alpha); s1 = sind(alpha);
    c2 = cosd(beta); s2 = sind(beta);
    c3 = cosd(gamma); s3 = sind(gamma);
    r = [c1*c2*c3-s1*s3 -c3*s1-c1*c2*s3 c1*s2; c1*s3+c2*c3*s1 c1*c3-c2*s1*s3 s1*s2; -c3*s2 s2*s3 c2]; %Z3 then Y then Z1
    
    %% Matrices transformation
    % Piezo-strain transformation
    d_rot = transform(d, r);
    d_rot_full = [d_rot(1,:,1) d_rot(2,:,1) d_rot(3,:,1); d_rot(1,:,2) d_rot(2,:,2) d_rot(3,:,2); d_rot(1,:,3) d_rot(2,:,3) d_rot(3,:,3)];
    d_rot_reduced = [d_rot_full(:,1) d_rot_full(:,5) d_rot_full(:,9) d_rot_full(:,2) d_rot_full(:,6) d_rot_full(:,3)];
    
    % T_sigma calculation
    Ts1 = [r(1,1)^2 r(2,1)^2 r(3,1)^2 2*r(1,1)*r(2,1) 2*r(2,1)*r(3,1) 2*r(3,1)*r(1,1)];
    Ts2 = [r(1,2)^2 r(2,2)^2 r(3,2)^2 2*r(1,2)*r(2,2) 2*r(2,2)*r(3,2) 2*r(3,2)*r(1,2)];
    Ts3 = [r(1,3)^2 r(2,3)^2 r(3,3)^2 2*r(1,3)*r(2,3) 2*r(2,3)*r(3,3) 2*r(3,3)*r(1,3)];
    Ts4 = [r(1,1)*r(1,2) r(2,1)*r(2,2) r(3,1)*r(3,2) r(1,1)*r(2,2)+r(1,2)*r(2,1) r(2,1)*r(3,2)+r(3,1)*r(2,2) r(3,1)*r(1,2)+r(3,2)*r(1,1)];
    Ts5 = [r(1,2)*r(1,3) r(2,2)*r(2,3) r(3,3)*r(3,2) r(2,3)*r(1,2)+r(1,3)*r(2,2) r(2,2)*r(3,3)+r(3,2)*r(2,3) r(3,2)*r(1,3)+r(1,2)*r(3,3)];
    Ts6 = [r(1,1)*r(1,3) r(2,3)*r(2,1) r(3,3)*r(3,1) r(1,3)*r(2,1)+r(1,1)*r(2,3) r(2,3)*r(3,1)+r(2,1)*r(3,3) r(3,3)*r(1,1)+r(3,1)*r(1,3)];
    Ts = [Ts1;Ts2;Ts3;Ts4;Ts5;Ts6];
    
    % T_epsilon calculation
    Te1 = [r(1,1)^2 r(2,1)^2 r(3,1)^2 r(1,1)*r(2,1) r(2,1)*r(3,1) r(3,1)*r(1,1)];
    Te2 = [r(1,2)^2 r(2,2)^2 r(3,2)^2 r(1,2)*r(2,2) r(2,2)*r(3,2) r(3,2)*r(1,2)];
    Te3 = [r(1,3)^2 r(2,3)^2 r(3,3)^2 r(1,3)*r(2,3) r(2,3)*r(3,3) r(3,3)*r(1,3)];
    Te4 = [2*r(1,1)*r(1,2) 2*r(2,1)*r(2,2) 2*r(3,1)*r(3,2) r(1,1)*r(2,2)+r(1,2)*r(2,1) r(2,1)*r(3,2)+r(3,1)*r(2,2) r(3,1)*r(1,2)+r(3,2)*r(1,1)];
    Te5 = [2*r(1,2)*r(1,3) 2*r(2,2)*r(2,3) 2*r(3,3)*r(3,2) r(2,3)*r(1,2)+r(1,3)*r(2,2) r(2,2)*r(3,3)+r(3,2)*r(2,3) r(3,2)*r(1,3)+r(1,2)*r(3,3)];
    Te6 = [2*r(1,1)*r(1,3) 2*r(2,3)*r(2,1) 2*r(3,3)*r(3,1) r(1,3)*r(2,1)+r(1,1)*r(2,3) r(2,3)*r(3,1)+r(2,1)*r(3,3) r(3,3)*r(1,1)+r(3,1)*r(1,3)];
    Te = [Te1;Te2;Te3;Te4;Te5;Te6];
    
    % Evaluation of rotated stiffness and rotated stress-charge [e] matrix
    c_rot = Ts^-1*c*Te;
    e_rot_ANSYS = transpose(d_rot_reduced*c_rot);
    d_rot_ANSYS = transpose(d_rot_reduced);
    d_rot_ANSYS = [d_rot_ANSYS(1,:);d_rot_ANSYS(2,:);d_rot_ANSYS(3,:);d_rot_ANSYS(5,:);d_rot_ANSYS(4,:);d_rot_ANSYS(6,:)];
    perm = [perm_x; perm_y; perm_z];
    perm_rot = abs(r*perm);
    
    %% Export to APDL commands
    % Permittivity
    p_join_4 = ['MP,PERX,',num2str(n+length(Cube)^3) ',', num2str(perm_rot(1))];
    p_join_5 = ['MP,PERY,',num2str(n+length(Cube)^3) ',', num2str(perm_rot(2))];
    p_join_6 = ['MP,PERZ,',num2str(n+length(Cube)^3) ',', num2str(perm_rot(3))];
    
    element_id_c_poled = ['TB,ANEL,',num2str(n+length(Cube)^3),',1,,0'];
    c_join1 = ['TBDATA,1,',num2str(c_rot(1,1)) ',',num2str(c_rot(1,2)) ',',num2str(c_rot(1,3))];
    c_join2 = ['TBDATA,4,',num2str(c_rot(1,4)) ',',num2str(c_rot(1,5)) ',',num2str(c_rot(1,6))];
    c_join3 = ['TBDATA,7,',num2str(c_rot(2,2)) ',',num2str(c_rot(2,3)) ',',num2str(c_rot(2,4))];
    c_join4 = ['TBDATA,10,',num2str(c_rot(2,5)) ',',num2str(c_rot(2,6)) ',',num2str(c_rot(3,3))];
    c_join5 = ['TBDATA,13,',num2str(c_rot(3,4)) ',',num2str(c_rot(3,5)) ',',num2str(c_rot(3,6))];
    c_join6 = ['TBDATA,16,',num2str(c_rot(4,4)) ',',num2str(c_rot(4,5)) ',',num2str(c_rot(4,6))];
    c_join7 = ['TBDATA,19,',num2str(c_rot(5,5)) ',',num2str(c_rot(5,6)) ',',num2str(c_rot(6,6))];
    
    element_id_e_poled = ['TB,PIEZ,',num2str(n+length(Cube)^3)];
    e_join1 = ['TBDATA,1,',num2str(e_rot_ANSYS(1,1)) ',',num2str(e_rot_ANSYS(1,2)) ',',num2str(e_rot_ANSYS(1,3))];
    e_join2 = ['TBDATA,4,',num2str(e_rot_ANSYS(2,1)) ',',num2str(e_rot_ANSYS(2,2)) ',',num2str(e_rot_ANSYS(2,3))];
    e_join3 = ['TBDATA,7,',num2str(e_rot_ANSYS(3,1)) ',',num2str(e_rot_ANSYS(3,2)) ',',num2str(e_rot_ANSYS(3,3))];
    e_join4 = ['TBDATA,10,',num2str(e_rot_ANSYS(4,1)) ',',num2str(e_rot_ANSYS(4,2)) ',',num2str(e_rot_ANSYS(4,3))];
    e_join5 = ['TBDATA,13,',num2str(e_rot_ANSYS(5,1)) ',',num2str(e_rot_ANSYS(5,2)) ',',num2str(e_rot_ANSYS(5,3))];
    e_join6 = ['TBDATA,16,',num2str(e_rot_ANSYS(6,1)) ',',num2str(e_rot_ANSYS(6,2)) ',',num2str(e_rot_ANSYS(6,3))];
    
    fprintf(fid, '%s', p_join_4);
    fprintf(fid, '\n');
    fprintf(fid, '%s', p_join_5);
    fprintf(fid, '\n');
    fprintf(fid, '%s', p_join_6);
    fprintf(fid, '\n');
    fprintf(fid, '%s', element_id_c_poled);
    fprintf(fid, '\n');
    fprintf(fid, c_join1);
    fprintf(fid, '\n');
    fprintf(fid, c_join2);
    fprintf(fid, '\n');
    fprintf(fid, c_join3);
    fprintf(fid, '\n');
    fprintf(fid, c_join4);
    fprintf(fid, '\n');
    fprintf(fid, c_join5);
    fprintf(fid, '\n');
    fprintf(fid, c_join6);
    fprintf(fid, '\n');
    fprintf(fid, c_join7);
    fprintf(fid, '\n');
    fprintf(fid, '%s', element_id_e_poled);
    fprintf(fid, '\n');
    fprintf(fid, e_join1);
    fprintf(fid, '\n');
    fprintf(fid, e_join2);
    fprintf(fid, '\n');
    fprintf(fid, e_join3);
    fprintf(fid, '\n');
    fprintf(fid, e_join4);
    fprintf(fid, '\n');
    fprintf(fid, e_join5);
    fprintf(fid, '\n');
    fprintf(fid, e_join6);
    fprintf(fid, '\n');

end

fclose('all')
toc




     




