%MeshGen.m
%=========================================================================%
% MatLab script to generate LS-DYNA mesh file for a radially symetric mesh
% for explosive (part 1) and air (part 2). Written by Samuel E. Rigby
%=========================================================================%
clear
clc

% INPUTS
% User defined mesh controls
charge = input('Enter charge radius (m): ');
airdomain = input('Enter domain length [includes charge radius] (m): ');

radial = input('Enter number of elements along radius of air domain: ');
A = 0;
while A~=2
    circumferential = input('Enter number of circumferential elements (must be even): ');
    fac = factor(circumferential);
    A = fac(1);
    if A~=2
        disp('Circumferential elements must be even')
        circumferential = input('Enter number of circumferential elements (must be even): ');
        fac = factor(circumferential);
        A = fac(1);
    end
end
clearvars A

innerchargediag = charge / 2;
innercharge = sqrt((innerchargediag^2)/2);
%hypotenuse of square mesh of charge = half radius of charge

innerchargeelms = circumferential / 2; %elements in innercharge
radialnodes = circumferential+1;
innerchargenodes = innerchargeelms+1;
outerchargenodes = round(((charge-innercharge)/innercharge)*innerchargenodes);
chargenodes = innerchargenodes + outerchargenodes;

xos = innercharge*ones(radialnodes,1);
xos(1:innerchargenodes) = linspace(0,innercharge,innerchargenodes);
yos = flipud(xos);

theta = linspace(0,pi/2,radialnodes);
theta = theta';
xends = charge*(cos((pi/2)-theta));
yends = flipud(xends);

for i = 1:size(xos);
    diagx(i,1:outerchargenodes) = linspace(xos(i),xends(i),outerchargenodes);
    diagy(i,1:outerchargenodes) = linspace(yos(i),yends(i),outerchargenodes);
end

xxos = linspace(0,innercharge,innerchargenodes);
yyos = linspace(0,innercharge,innerchargenodes);
[gridx,gridy] = meshgrid(xxos,yyos);
gridx(:,end) = [];
gridy(:,end) = [];
gridx(end,:) = [];
gridy(end,:) = [];

diagX = reshape(diagx,numel(diagx),1);
diagY = reshape(diagy,numel(diagy),1);

gridX = reshape(gridx,numel(gridx),1);
gridY = reshape(gridy,numel(gridy),1);

numsgrid = numel(gridX);
numsdiag = numel(diagX);
nodecount = numsgrid+numsdiag;
totalchargeelms = ((outerchargenodes-1)*circumferential) + innerchargeelms^2;

X = zeros(numsgrid+numsdiag,1);
Y = X;
X(1:numsgrid) = gridX;
X(numsgrid+1:end) = diagX;
Y(1:numsgrid) = gridY;
Y(numsgrid+1:end) = diagY;

str = [num2str(circumferential) 'circumferential' num2str(radial) 'radial_mesh.k']; 
filename = str;
fp = fopen(filename,'w');

fprintf(fp,'*KEYWORD\n');
fprintf(fp,'*NODE\n');
fprintf(fp,'$#   nid               x               y               z      tc      rc\n');
% NODES___________________________________________________________
for i = 1:nodecount
    nodenum = i;
    xloc = X(i);
    yloc = Y(i);
    fprintf(fp,'%8s',num2str(nodenum));
    fprintf(fp,'%16s%16s',num2str(xloc),num2str(yloc));
    fprintf(fp,'           0.000       0       0\n');
end

% AIR MESH
l = logspace(log10(charge),log10(airdomain),radial+1);
% l = linspace(charge,airdomain,airdomainelms+1);
l(1) = [];
theta = linspace(pi/2,0,radialnodes);

[R,THETA] = meshgrid(l,theta);
[Xair,Yair] = pol2cart(THETA,R);

XX = reshape(Xair,numel(Xair),1);
YY = reshape(Yair,numel(Yair),1);

for i = 1:numel(Xair)
    nodenum = nodecount + i;
    xloc = XX(i);
    yloc = YY(i);
    fprintf(fp,'%8s',num2str(nodenum));
    fprintf(fp,'%16s%16s',num2str(xloc),num2str(yloc));
    fprintf(fp,'           0.000       0       0\n');
end

% ELEMENTS________________________________________________________
fprintf(fp,'*ELEMENT_SHELL\n');
fprintf(fp,'$#   eid     pid      n1      n2      n3      n4      n5      n6      n7      n8\n');

gridelecount = innerchargeelms-1;
gridnodecount = innerchargenodes-1;

% ELEMENTS IN MAIN GRID OF EXPLOSIVE
for i = 1:gridelecount;
    for j = 1:gridelecount;
        elenum = ((gridelecount * i) + j - gridelecount);
        n1 = ((i - 1) * gridnodecount) + j + 2 + gridelecount;
        n2 = ((i - 1) * gridnodecount) + j + 1;
        n3 = ((i - 1) * gridnodecount) + j;
        n4 = ((i - 1) * gridnodecount) + j + 1 + gridelecount;
        fprintf(fp,'%8s%8s',num2str(elenum),num2str(1));
        fprintf(fp,'%8s%8s%8s%8s',num2str(n1),num2str(n2),num2str(n3),num2str(n4));
        fprintf(fp,'       0       0       0       0\n');
    end
end

% ELEMENTS ALONG TOP ROW OF EXPLOSIVE GRID
newelstart = (gridelecount^2)+1;
newelend = (innerchargeelms^2)-(innerchargeelms-1);
diff = newelend - newelstart;
for i = 1:diff+1
    elenum = newelstart + (i-1);
    n1 = (gridnodecount^2) + (i+1);
    if i == diff+1
    n2 = n1 + 1;
    else
    n2 = gridnodecount * (i+1);
    end
    n3 = gridnodecount * i;
    n4 = (gridnodecount^2) + i;
    fprintf(fp,'%8s%8s',num2str(elenum),num2str(1));
    fprintf(fp,'%8s%8s%8s%8s',num2str(n1),num2str(n2),num2str(n3),num2str(n4));
    fprintf(fp,'       0       0       0       0\n');
end

% ELEMENTS ALONG RHS OF EXPLOSIVE GRID
newelstart = newelend + 1;
newelend = (innerchargeelms^2);
diff = newelend - newelstart;
for i = 1:diff+1
    elenum = newelstart + (i-1);
    n1 = (gridnodecount^2) + (innerchargenodes) + i;
    n2 = (gridnodecount^2) + (innerchargenodes) + i + 1;
    n3 = (gridnodecount^2) - i;
    n4 = (gridnodecount^2) - (i - 1);
    fprintf(fp,'%8s%8s',num2str(elenum),num2str(1));
    fprintf(fp,'%8s%8s%8s%8s',num2str(n1),num2str(n2),num2str(n3),num2str(n4));
    fprintf(fp,'       0       0       0       0\n');
end

% REMAINING ELEMENTS
for i = 1:outerchargenodes - 1;
    for j = 1:circumferential;
        elenum = (innerchargeelms^2) + ((i-1) * circumferential) + j;
        n1 = (innerchargenodes^2) + ((i-1) * radialnodes) + (j+1);
        n2 = (gridnodecount^2) + ((i-1) * radialnodes) + (j+1);
        n3 = n2 - 1;
        n4 = n1 - 1;
        fprintf(fp,'%8s%8s',num2str(elenum),num2str(1));
        fprintf(fp,'%8s%8s%8s%8s',num2str(n1),num2str(n2),num2str(n3),num2str(n4));
        fprintf(fp,'       0       0       0       0\n');
    end
end

elenum = totalchargeelms;
airnodes = (radial*circumferential);
for i = 1:radial;
    for j = 1:circumferential;
        elenum = elenum+1;
        n1 = (nodecount) + ((i-1) * radialnodes) + (j+1);
        n2 = (nodecount) + ((i-2) * radialnodes) + (j+1);
        n3 = n2 - 1;
        n4 = n1 - 1;
        fprintf(fp,'%8s%8s',num2str(elenum),num2str(2));
        fprintf(fp,'%8s%8s%8s%8s',num2str(n1),num2str(n2),num2str(n3),num2str(n4));
        fprintf(fp,'       0       0       0       0\n');
    end
end

% Y BOUNDARY NODES
boundnodeschargeY = find(Y<1E-15);
boundnodesairY = find(Yair<1E-15);
boundnodesairY = boundnodesairY + numel(Y);
totalboundaryY = numel(boundnodeschargeY) + numel(boundnodesairY);

rows = ceil(totalboundaryY / 8);
numboundary = rows * 8;
boundarynodes = zeros(numboundary,1);

%Set up boundary nodes as a vector
boundarynodes(1:numel(boundnodeschargeY)) = boundnodeschargeY;
boundarynodes(numel(boundnodeschargeY)+1:totalboundaryY) = boundnodesairY;

%Put vector into nx8 matrix
bound = zeros(rows,8);

fprintf(fp,'*BOUNDARY_SPC_SET\n');
fprintf(fp,'$#    nsid       cid      dofx      dofy      dofz     dofrx     dofry     dofrz\n');
fprintf(fp,'         1         0         0         1         1         1         1         1\n');
fprintf(fp,'*SET_NODE_LIST_TITLE\n');
fprintf(fp,'Boundary Nodes Along Bottom\n');
fprintf(fp,'$#     sid       da1       da2       da3       da4    solver\n');
fprintf(fp,'         1     0.000     0.000     0.000     0.000MECH\n');
fprintf(fp,'$#    nid1      nid2      nid3      nid4      nid5      nid6      nid7      nid8\n');

for i = 1:rows
    bound(i,:) = boundarynodes((8*(i-1))+1:8*i);
    fprintf(fp,'%10s%10s%10s%10s%10s%10s%10s%10s\n',num2str(bound(i,1)),...
        num2str(bound(i,2)),num2str(bound(i,3)),num2str(bound(i,4)),...
        num2str(bound(i,5)),num2str(bound(i,6)),num2str(bound(i,7)),...
        num2str(bound(i,8)));
end

% X BOUNDARY NODES
boundnodeschargeX = find(X<1E-15);
boundnodesairX = find(Xair<1E-15);
boundnodesairX = boundnodesairX + numel(X);
totalboundaryX = numel(boundnodeschargeX) + numel(boundnodesairX);

rowsX = ceil(totalboundaryX / 8);
numboundaryX = rowsX * 8;
boundarynodesX = zeros(numboundaryX,1);

%Set up boundary nodes as a vector
boundarynodesX(1:numel(boundnodeschargeX)) = boundnodeschargeX;
boundarynodesX(numel(boundnodeschargeX)+1:totalboundaryX) = boundnodesairX;

%Put vector into nx8 matrix
boundX = zeros(rows,8);

fprintf(fp,'*BOUNDARY_SPC_SET\n');
fprintf(fp,'$#    nsid       cid      dofx      dofy      dofz     dofrx     dofry     dofrz\n');
fprintf(fp,'         2         0         1         0         1         1         1         1\n');
fprintf(fp,'*SET_NODE_LIST_TITLE\n');
fprintf(fp,'Boundary Nodes Along Symm Axis\n');
fprintf(fp,'$#     sid       da1       da2       da3       da4    solver\n');
fprintf(fp,'         2     0.000     0.000     0.000     0.000MECH\n');
fprintf(fp,'$#    nid1      nid2      nid3      nid4      nid5      nid6      nid7      nid8\n');

for i = 1:rows
    boundX(i,:) = boundarynodesX((8*(i-1))+1:8*i);
    fprintf(fp,'%10s%10s%10s%10s%10s%10s%10s%10s\n',num2str(boundX(i,1)),...
        num2str(boundX(i,2)),num2str(boundX(i,3)),num2str(boundX(i,4)),...
        num2str(boundX(i,5)),num2str(boundX(i,6)),num2str(boundX(i,7)),...
        num2str(boundX(i,8)));
end
fprintf(fp,'*END');

nodeadd = size(X,1);
aspectlength = ( (XX(n1-nodeadd) - XX(n2-nodeadd))^2 +...
    (YY(n1-nodeadd) - YY(n2-nodeadd))^2 )^0.5;
aspectheight = ( (YY(n4-nodeadd) - YY(n1-nodeadd))^2 +...
    (XX(n4-nodeadd) - XX(n1-nodeadd))^2 )^0.5;
clc
disp('Meshing complete')
aspectratio = min(aspectlength,aspectheight)/max(aspectlength,aspectheight);
if aspectratio < 0.1
    disp('WARNING: Aspect ratio of air mesh < 0.1')
    disp('Please increase circumferential elements for stable calculation')
end