% // We need this to start
import com.comsol.model.util.*
% // You should be in the directory where this file is located;
% // and/or change the location and file name. 
%model = mphload('Hip');
%model = mphload('small_tet');
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

% //Based on the stats info, we see we have 4 element types
% //this gives element connectivity for the 4 element types; note the transpose
idx_tetNodeConnectivity = info.elements.tet.nodes(:,:)';
idx_edgNodeConnectivity = info.elements.edg.nodes(:,:)';
idx_triNodeConnectivity = info.elements.tri.nodes(:,:)';
idx_vtxNodeConnectivity = info.elements.vtx.nodes(:,:)';

% These are all a 0-based connectivity info. Convert to a 1=based
idx_tetNodeConnectivity = idx_tetNodeConnectivity + 1
idx_edgNodeConnectivity = idx_edgNodeConnectivity + 1
idx_triNodeConnectivity = idx_triNodeConnectivity + 1
idx_vtxNodeConnectivity = idx_vtxNodeConnectivity + 1

% TETS
numTetElements = size(idx_tetNodeConnectivity, 1);
tetElementConnectivity = zeros(numTetElements, 11);
for i = 1:numTetElements
	tetElementConnectivity(i, 1) =  i;
	tetElementConnectivity(i, 2) =  idx_tetNodeConnectivity(i, 1);
	tetElementConnectivity(i, 3) =  idx_tetNodeConnectivity(i, 2);
	tetElementConnectivity(i, 4) =  idx_tetNodeConnectivity(i, 3);
	tetElementConnectivity(i, 5) =  idx_tetNodeConnectivity(i, 4);
	tetElementConnectivity(i, 6) =  idx_tetNodeConnectivity(i, 5);
	tetElementConnectivity(i, 7) =  idx_tetNodeConnectivity(i, 6);
	tetElementConnectivity(i, 8) =  idx_tetNodeConnectivity(i, 7);
	tetElementConnectivity(i, 9) =  idx_tetNodeConnectivity(i, 8);
	tetElementConnectivity(i, 10) = idx_tetNodeConnectivity(i, 9);
	tetElementConnectivity(i, 11) = idx_tetNodeConnectivity(i, 10);
end
writematrix(tetElementConnectivity, 'ElementConnectivityTET','Delimiter',',')
type ElementConnectivityTET.txt

% 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



% TRI
numTriElements = size(idx_triNodeConnectivity, 1);
triElementConnectivity = zeros(numTriElements, 7);
for i = 1:numTriElements
	triElementConnectivity(i, 1) =  i;
	triElementConnectivity(i, 2) =  idx_triNodeConnectivity(i, 1);
	triElementConnectivity(i, 3) =  idx_triNodeConnectivity(i, 2);
	triElementConnectivity(i, 4) =  idx_triNodeConnectivity(i, 3);
	triElementConnectivity(i, 5) =  idx_triNodeConnectivity(i, 4);
	triElementConnectivity(i, 6) =  idx_triNodeConnectivity(i, 5);
	triElementConnectivity(i, 7) =  idx_triNodeConnectivity(i, 6);
end
writematrix(triElementConnectivity, 'ElementConnectivityTRI','Delimiter', ',')
type ElementConnectivityTRI.txt

% EDG
numEdgElements = size(idx_edgNodeConnectivity, 1);
edgElementConnectivity = zeros(numEdgElements, 4);
for i = 1:numEdgElements
	edgElementConnectivity(i, 1) =  i;
	edgElementConnectivity(i, 2) =  idx_edgNodeConnectivity(i, 1);
	edgElementConnectivity(i, 3) =  idx_edgNodeConnectivity(i, 2);
	edgElementConnectivity(i, 4) =  idx_edgNodeConnectivity(i, 3);
end
writematrix(edgElementConnectivity, 'ElementConnectivityEDG','Delimiter', ',')
type ElementConnectivityEDG.txt

% VTX
numVtxElements = size(idx_vtxNodeConnectivity, 1);
vtxElementConnectivity = zeros(numVtxElements, 2);
for i = 1:numVtxElements
	vtxElementConnectivity(i, 1) =  i;
	vtxElementConnectivity(i, 2) =  idx_vtxNodeConnectivity(i, 1);
end
writematrix(vtxElementConnectivity, 'ElementConnectivityVTX','Delimiter', ',')
type ElementConnectivityVTX.txt

% This is everthing we need for the input file...
% This is what we need for the results

% DISPLACEMENTS

DispData = mpheval(model, {'u', 'v', 'w' }, 'refine', 2);

numDisplacements = size(DispData.d1', 1);
Disp = [DispData.d1', DispData.d2', DispData.d3']
Displacements = zeros(numDisplacements, 4);
for i = 1:numDisplacements;
    Displacements(i, 1) = i;
    Displacements(i, 2) = Disp(i, 1);
    Displacements(i, 3) = Disp(i, 2);
    Displacements(i, 4) = Disp(i, 3);
end
writematrix(Displacements, 'Displacements','Delimiter', ',')
type Displacements.txt


% STRESS

StressData = mpheval(model, {'solid.sGpxx', 'solid.sGpyy', 'solid.sGpzz', 'solid.sGpxy', 'solid.sGpxz', 'solid.sGpyz'}, 'refine', 2)

numStress = size(StressData.d1', 1);
Stresses = [StressData.d1', StressData.d2', StressData.d3', StressData.d4', StressData.d5', StressData.d6']
Stress = zeros(numStress, 7);
for i = 1:numStress
    Stress(i, 1) = i;
    Stress(i, 2) = Stresses(i, 1);
    Stress(i, 3) = Stresses(i, 2);
    Stress(i, 4) = Stresses(i, 3);
    Stress(i, 5) = Stresses(i, 4);
    Stress(i, 6) = Stresses(i, 5);
    Stress(i, 7) = Stresses(i, 6);
end
writematrix(Stress, 'Stress','Delimiter', ',')
type Stress.txt

% STRAIN

StrainData = mpheval(model, {'solid.eeGpXX', 'solid.eeGpYY', 'solid.eeGpZZ', 'solid.eeGpXY', 'solid.eeGpXY', 'solid.eeGpYZ'}, 'refine', 2)

numStrain = size(StrainData.d1', 1);
Strains = [StrainData.d1', StrainData.d2', StrainData.d3', StrainData.d4', StrainData.d5', StrainData.d6']
Strain = zeros(numStrain, 7);
for i = 1:numStrain
    Strain(i, 1) = i;
    Strain(i, 2) = Strains(i, 1);
    Strain(i, 3) = Strains(i, 2);
    Strain(i, 4) = Strains(i, 3);
    Strain(i, 5) = Strains(i, 4);
    Strain(i, 6) = Strains(i, 5);
    Strain(i, 7) = Strains(i, 6);
end
writematrix(Strain, 'Strain','Delimiter', ',')
type Strain.txt