% utmc.m
% Uniform tetrahedral meshes for cubes (3-dim rectangles)
% Saving info to disk as soon as generated
% Jiangguo (James) Liu, ColoState;  Haiying Wang, MTU;  10/2009--04/2013

function  [] = utmc()
  lblCrnr = zeros(1,4);
  lblLast = zeros(1,6);
  lblThis = zeros(1,6);
  lblVrtx = zeros(1,8);

  disp('X-direction: ');
  xa = input('  Left endpoint  (xa) = ');
  xb = input('  Right endpoint (xb) = ');
  nx = input('  Number of partitions (nx) = ');

  disp('Y-direction: ');
  yc = input('  Left endpoint  (yc) = ');
  yd = input('  Right endpoint (yd) = ');
  ny = input('  Number of partitions (ny) = ');

  disp('Z-direction: ');
  ze = input('  Left endpoint  (ze) = ');
  zf = input('  Right endpoint (zf) = ');
  nz = input('  Number of partitions (nz) = ');

  x = linspace(xa,xb,nx+1);
  hx = (xb-xa)/nx;

  y = linspace(yc,yd,ny+1);
  hy = (yd-yc)/ny;

  z = linspace(ze,zf,nz+1);
  hz = (zf-ze)/nz;

  % Saving info to disk as soon as generated
  fp = fopen('mTetraMesh.dat', 'wt');
  fp0 = fopen('mTetraMesh.node', 'wt');
  fp2 = fopen('mTetraMesh.face', 'wt');
  fp3 = fopen('mTetraMesh.ele', 'wt');

  % Saving the numbers of nodes/faces/elements/boundarypieces
  % By default, all strating labels are 1
  numNds = (nx+1)*(ny+1)*(nz+1);
  numFcs = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1) + 2*ny*nz*(nx+1) + 6*nx*ny*nz;
  numElts = 6*nx*ny*nz;
  numPcsBndry = 6;
  fprintf(fp, '%7d  %7d  %7d  %7d\n', numNds, numFcs, numElts, numPcsBndry);
  fprintf(fp, '%2d  %2d  %2d  %2d\n', 1, 1, 1, 1);
  fprintf(fp0, '%7d  %2d  %2d  %2d\n', numNds, 3, 0, 0);
  fprintf(fp2, '%7d  %2d\n', numFcs, 1);
  fprintf(fp3, '%7d  %2d  %2d\n', numElts, 4, 0);

  % Saving all node coordinates
  labeld = 1;
  for k = 1 : nz+1 
    for j = 1 : ny+1 
      for i = 1 : nx+1 
        if (i==1 | i==nx+1 | j==1 | j==ny+1 | k==1 | k==nz+1) 
          isBndryNd = 1;
        else 
          isBndryNd = 0;
        end
        fprintf(fp, '%7d   %12.6f  %12.6f  %12.6f  %2d\n', labeld, ...
          x(i), y(j), z(k), isBndryNd);
        fprintf(fp0, '%7d   %12.6f  %12.6f  %12.6f\n', labeld, ...
          x(i), y(j), z(k));
        labeld = labeld + 1;
      end
    end
  end
  clear  x  y  z;  % Not really needed

  % Generating and saving info on the block faces:  XY-directions
  labelf = 0;
  for k = 1 : nz+1 
    isBndryFc = 0;
    if (k==1)  isBndryFc = 1;  end
    if (k==nz+1)  isBndryFc = 2;  end
    for j = 1 : ny 
      for i = 1 : nx 
        lblLast = cmptGlbLblsTetraBlk(i, j, k-1, nx, ny, nz);
        lblThis = cmptGlbLblsTetraBlk(i, j, k, nx, ny, nz);

        lblCrnr(1) = glbLblNd(i, j, k, nx, ny);
        lblCrnr(2) = glbLblNd(i+1, j, k, nx, ny);
        lblCrnr(3) = glbLblNd(i+1, j+1, k, nx, ny);
        lblCrnr(4) = glbLblNd(i, j+1, k, nx, ny);

        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+1, ...
          lblCrnr(1), lblCrnr(3), lblCrnr(4), isBndryFc, ...
          lblLast(3), lblThis(1));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+2, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(3), isBndryFc, ...
          lblLast(6), lblThis(4));

        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+1, ...
          lblCrnr(1), lblCrnr(3), lblCrnr(4), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+2, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(3), isBndryFc);

        labelf = labelf + 2;
      end
    end
  end

  % Generating and saving info on the block faces:  XZ-directions
  labelf = 2*nx*ny*(nz+1);
  for j = 1 : ny+1 
    isBndryFc = 0;
    if (j==1)  isBndryFc = 3;  end
    if (j==ny+1)  isBndryFc = 4;  end
    for k = 1 : nz 
      for i = 1 : nx 
        lblLast = cmptGlbLblsTetraBlk(i, j-1, k, nx, ny, nz);
        lblThis = cmptGlbLblsTetraBlk(i, j, k, nx, ny, nz);

        lblCrnr(1) = glbLblNd(i, j, k, nx, ny);
        lblCrnr(2) = glbLblNd(i+1, j, k, nx, ny);
        lblCrnr(3) = glbLblNd(i+1, j, k+1, nx, ny);
        lblCrnr(4) = glbLblNd(i, j, k+1, nx, ny);

        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+1, ... 
          lblCrnr(1), lblCrnr(3), lblCrnr(4), isBndryFc, ...
          lblLast(3), lblThis(5));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+2, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(3), isBndryFc, ...
          lblLast(2), lblThis(4));

        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+1, ... 
          lblCrnr(1), lblCrnr(3), lblCrnr(4), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+2, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(3), isBndryFc);

        labelf = labelf + 2;
      end
    end
  end

  % Generating and saving info on the block faces:  YZ-directions
  labelf = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1);
  for i = 1 : nx+1 
    isBndryFc = 0;
    if (i==1)  isBndryFc = 5;  end
    if (i==nx+1)  isBndryFc = 6;  end
    for k = 1 : nz 
      for j = 1 : ny 
        lblLast = cmptGlbLblsTetraBlk(i-1, j, k, nx, ny, nz);
        lblThis = cmptGlbLblsTetraBlk(i, j, k, nx, ny, nz);

        lblCrnr(1) = glbLblNd(i, j, k, nx, ny);
        lblCrnr(2) = glbLblNd(i, j+1, k, nx, ny); 
        lblCrnr(3) = glbLblNd(i, j+1, k+1, nx, ny);
        lblCrnr(4) = glbLblNd(i, j, k+1, nx, ny);

        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+1, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(4), isBndryFc, ...
          lblLast(4), lblThis(1));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+2, ...
          lblCrnr(2), lblCrnr(3), lblCrnr(4), isBndryFc, ...
          lblLast(6), lblThis(3));

        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+1, ...
          lblCrnr(1), lblCrnr(2), lblCrnr(4), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+2, ...
          lblCrnr(2), lblCrnr(3), lblCrnr(4), isBndryFc);

        labelf = labelf + 2;
      end
    end
  end

  % Generating and saving info on the interior faces
  labelf = 2*nx*ny*(nz+1) + 2*nx*nz*(ny+1) + 2*ny*nz*(nx+1);
  isBndryFc = 0;
  for k = 1 : nz 
    for j = 1 : ny 
      for i = 1 : nx 
        lblVrtx = cmptGlbLblsVrtxBlk(i, j, k, nx, ny);
        lblThis = cmptGlbLblsTetraBlk(i, j, k, nx, ny, nz);
        
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+1, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(5), isBndryFc, ...
          lblThis(1), lblThis(5));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+2, ...
          lblVrtx(3), lblVrtx(5), lblVrtx(7), isBndryFc, ...
          lblThis(2), lblThis(6));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+3, ...
          lblVrtx(3), lblVrtx(4), lblVrtx(5), isBndryFc, ... 
          lblThis(1), lblThis(2));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+4, ...
          lblVrtx(4), lblVrtx(5), lblVrtx(7), isBndryFc, ...
          lblThis(2), lblThis(3));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+5, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(6), isBndryFc, ...
          lblThis(4), lblThis(5));
        fprintf(fp, '%7d   %7d  %7d  %7d  %2d  %7d  %7d\n', labelf+6, ... 
          lblVrtx(3), lblVrtx(5), lblVrtx(6), isBndryFc, ...
          lblThis(5), lblThis(6));

        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+1, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(5), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+2, ...
          lblVrtx(3), lblVrtx(5), lblVrtx(7), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+3, ...
          lblVrtx(3), lblVrtx(4), lblVrtx(5), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+4, ... 
          lblVrtx(4), lblVrtx(5), lblVrtx(7), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+5, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(6), isBndryFc);
        fprintf(fp2, '%7d   %7d  %7d  %7d  %2d\n', labelf+6, ...
          lblVrtx(3), lblVrtx(5), lblVrtx(6), isBndryFc);

        labelf = labelf + 6;
      end
    end
  end

  % Generating and saving element info
  labele = 1;
  for k = 1 : nz 
    for j = 1 : ny 
      for i = 1 : nx 
        lblVrtx = cmptGlbLblsVrtxBlk(i, j, k, nx, ny);

        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(4), lblVrtx(5));
        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele+1, ...
          lblVrtx(3), lblVrtx(4), lblVrtx(5), lblVrtx(7));
        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele+2, ...
          lblVrtx(4), lblVrtx(5), lblVrtx(7), lblVrtx(8));
        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele+3, ...
          lblVrtx(1), lblVrtx(2), lblVrtx(3), lblVrtx(6));
        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele+4, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(5), lblVrtx(6));
        fprintf(fp, '%7d   %7d  %7d  %7d  %7d\n', labele+5, ...
          lblVrtx(3), lblVrtx(5), lblVrtx(6), lblVrtx(7));

        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele, ... 
          lblVrtx(1), lblVrtx(3), lblVrtx(4), lblVrtx(5));
        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele+1, ...
          lblVrtx(3), lblVrtx(4), lblVrtx(5), lblVrtx(7));
        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele+2, ... 
          lblVrtx(4), lblVrtx(5), lblVrtx(7), lblVrtx(8));
        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele+3, ...
          lblVrtx(1), lblVrtx(2), lblVrtx(3), lblVrtx(6));
        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele+4, ...
          lblVrtx(1), lblVrtx(3), lblVrtx(5), lblVrtx(6));
        fprintf(fp3, '%7d   %7d  %7d  %7d  %7d\n', labele+5, ...
          lblVrtx(3), lblVrtx(5), lblVrtx(6), lblVrtx(7));

        labele = labele + 6;
      end
    end
  end

  % Generating and saving boundary info
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 1, 0.0, 0.0, -1.0);
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 2, 0.0, 0.0, 1.0);
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 3, 0.0, -1.0, 0.0);
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 4, 0.0, 1.0, 0.0);
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 5, -1.0, 0.0, 0.0);
  fprintf(fp, '%2d  %7.4f  %7.4f  %7.4f\n', 6, 1.0, 0.0, 0.0);

  fclose(fp);
  fclose(fp0);
  fclose(fp2);
  fclose(fp3);
  disp('Mesh data successfully written to files mTetraMesh.****');
  disp('Press any key to terminate...');
  pause;
return;


% Subfunction: 
% Computing the global labels of the 6 tetrahedra inside a given block
function  lblTetra = cmptGlbLblsTetraBlk(i, j, k, nx, ny, nz) 
  lblTetra = zeros(1,6);
  if (i<1 | i>nx | j<1 | j>ny | k<1 | k>nz) 
    for m = 1 : 6 
      lblTetra(m) = -1;
    end
  else
    n = ((k-1)*ny*nx+(j-1)*nx+i-1)*6;
    for m = 1 : 6 
      lblTetra(m) = m+n;
    end
  end
return;


% Subfunction: 
% Computing the global labels of the 8 vertices of a given block
function  lblVrtx = cmptGlbLblsVrtxBlk(i, j, k, nx, ny) 
  lblVrtx = zeros(1,8);
  lblVrtx(1) = (k-1)*(ny+1)*(nx+1) + (j-1)*(nx+1) + i;
  lblVrtx(2) = lblVrtx(1) + 1;
  lblVrtx(3) = lblVrtx(1) + nx+2;
  lblVrtx(4) = lblVrtx(1) + nx+1;
  lblVrtx(5) = lblVrtx(1) + (ny+1)*(nx+1);
  lblVrtx(6) = lblVrtx(5) + 1;
  lblVrtx(7) = lblVrtx(5) + nx+2;
  lblVrtx(8) = lblVrtx(5) + nx+1;
return;


% Subfunction: 
% The linear global label of a node 
function  lblNd = glbLblNd(i, j, k, nx, ny) 
  lblNd = (k-1)*(nx+1)*(ny+1) + (j-1)*(nx+1) + i;
return;

