% // We need this to start
import com.comsol.model.util.*
% // You should be in the directory where the .mph file is located;
% // and/or change the location and file name. 

%model = mphload('OneTet');

% // This launches the Comsol app with this file
mphlaunch;

% // This runs the file
mphrun(model);

% // MESH INFO: Nodes

stats = mphmeshstats(model)
[stats, data] = mphmeshstats(model);
info = mphxmeshinfo(model)

%// Get the number of degrees of freedom and the node coordinates
dofs = info.ndofs
dof_coords = info.dofs.coords';
coordNumber = size(dof_coords, 1);
Node_coords = info.nodes.coords';       % Transpose original array
nodeNumber = size(Node_coords, 1);      % Number of rows (points)

% Preallocate Nodes based on the number of points and the expected number of coordinates
Nodes = zeros(nodeNumber, 4);

% Check the number of columns in Newcoords
numCols = size(Node_coords, 2);
disp(['Size of Node_coords: ', num2str(nodeNumber), ' rows, ', num2str(numCols), ' columns']);

for i = 1:nodeNumber
    Nodes(i, 1) = i;  % Counter
    Nodes(i, 2) = Node_coords(i, min(1, numCols));  % First coordinate
    Nodes(i, 3) = Node_coords(i, min(2, numCols));  % Second coordinate
    Nodes(i, 4) = Node_coords(i, min(3, numCols));  % Third coordinate
end

writematrix(Nodes, 'Nodes','Delimiter',',')
type Nodes.txt
% This makes the node points 1-based, since we number them ourselves

% MESH INFO: Elements

% //We are looking only for the tet elements note the transpose
idx_tetNodeConnectivity = info.elements.tet.nodes(:,:)';

% Thhis is a 0-based connectivity info. Convert to a 1=based
idx_tetNodeConnectivity = idx_tetNodeConnectivity + 1

% ITS TET Connectivity
numTetElements = size(idx_tetNodeConnectivity, 1);
ITS_tetElementConnectivity = zeros(numTetElements, 11);
for i = 1:numTetElements
	ITS_tetElementConnectivity(i, 1) =  i;
	ITS_tetElementConnectivity(i, 2) =  idx_tetNodeConnectivity(i, 1);
	ITS_tetElementConnectivity(i, 3) =  idx_tetNodeConnectivity(i, 3);
	ITS_tetElementConnectivity(i, 4) =  idx_tetNodeConnectivity(i, 6);
	ITS_tetElementConnectivity(i, 5) =  idx_tetNodeConnectivity(i, 10);
	ITS_tetElementConnectivity(i, 6) =  idx_tetNodeConnectivity(i, 2);
	ITS_tetElementConnectivity(i, 7) =  idx_tetNodeConnectivity(i, 5);
	ITS_tetElementConnectivity(i, 8) =  idx_tetNodeConnectivity(i, 4);
	ITS_tetElementConnectivity(i, 9) =  idx_tetNodeConnectivity(i, 8);
	ITS_tetElementConnectivity(i, 10) = idx_tetNodeConnectivity(i, 9);
	ITS_tetElementConnectivity(i, 11) = idx_tetNodeConnectivity(i, 7);
end
writematrix(ITS_tetElementConnectivity, 'ITS_ElementConnectivityTET','Delimiter',',')
type ITS_ElementConnectivityTET.txt



% ---------------------------------------------
% RESULTS: Displacements

% Get node coordinates in the right format for mphinterp
node_coords_for_eval = Nodes(:, 2:4)'; % [3 x N] format

% Extract displacements at exact node locations
disp_u = mphinterp(model, 'u', 'coord', node_coords_for_eval);
disp_v = mphinterp(model, 'v', 'coord', node_coords_for_eval);
disp_w = mphinterp(model, 'w', 'coord', node_coords_for_eval);

% Assemble displacement data in ITS format
numDisplacements = size(Nodes, 1);
Disp = [disp_u', disp_v', disp_w'];
Displacements = zeros(numDisplacements, 4);
for i = 1:numDisplacements
    Displacements(i, 1) = i;
    Displacements(i, 2) = Disp(i, 1); % u displacement
    Displacements(i, 3) = Disp(i, 2); % v displacement
    Displacements(i, 4) = Disp(i, 3); % w displacement
end

% Save the displacement data
writematrix(Displacements, 'Displacements','Delimiter', ',');


% RESULTS: Stress

% Extract stress components at exact node locations
stress_sxx = mphinterp(model, 'solid.sxx', 'coord', node_coords_for_eval);
stress_syy = mphinterp(model, 'solid.syy', 'coord', node_coords_for_eval);
stress_szz = mphinterp(model, 'solid.szz', 'coord', node_coords_for_eval);
stress_sxy = mphinterp(model, 'solid.sxy', 'coord', node_coords_for_eval);
stress_sxz = mphinterp(model, 'solid.sxz', 'coord', node_coords_for_eval);
stress_syz = mphinterp(model, 'solid.syz', 'coord', node_coords_for_eval);

% Assemble stress data in ITS format
numStress = size(Nodes, 1);
Stresses = [stress_sxx', stress_syy', stress_szz', stress_sxy', stress_sxz', stress_syz'];
Stress = zeros(numStress, 7);
for i = 1:numStress
    Stress(i, 1) = i;
    Stress(i, 2) = Stresses(i, 1); % sxx
    Stress(i, 3) = Stresses(i, 2); % syy
    Stress(i, 4) = Stresses(i, 3); % szz
    Stress(i, 5) = Stresses(i, 4); % sxy
    Stress(i, 6) = Stresses(i, 5); % sxz
    Stress(i, 7) = Stresses(i, 6); % syz
end

% Save the stress data
writematrix(Stress, 'Stress','Delimiter', ',')


% RESULTS: Strain

% Extract strain components at exact node locations - BEST CHOICE
strain_exx = mphinterp(model, 'solid.eXX', 'coord', node_coords_for_eval);
strain_eyy = mphinterp(model, 'solid.eYY', 'coord', node_coords_for_eval);
strain_ezz = mphinterp(model, 'solid.eZZ', 'coord', node_coords_for_eval);
strain_exy = mphinterp(model, 'solid.eXY', 'coord', node_coords_for_eval);
strain_exz = mphinterp(model, 'solid.eXZ', 'coord', node_coords_for_eval);
strain_eyz = mphinterp(model, 'solid.eYZ', 'coord', node_coords_for_eval);

% Assemble strain data in ITS format (matching stress format)
numStrain = size(Nodes, 1);
Strains = [strain_exx', strain_eyy', strain_ezz', strain_exy', strain_exz', strain_eyz'];
Strain = zeros(numStrain, 7);
for i = 1:numStrain
    Strain(i, 1) = i;
    Strain(i, 2) = Strains(i, 1); % eXX (εxx)
    Strain(i, 3) = Strains(i, 2); % eYY (εyy)
    Strain(i, 4) = Strains(i, 3); % eZZ (εzz)
    Strain(i, 5) = Strains(i, 4); % eXY (εxy)
    Strain(i, 6) = Strains(i, 5); % eXZ (εxz)
    Strain(i, 7) = Strains(i, 6); % eYZ (εyz)
end

writematrix(Strain, 'Strain','Delimiter', ',');