// tri_split.java  split each triangle into four similar triangles
class tri_split
{
  // tri_split  split triangles, each one becomes four triangles
  //            similar to the original
  //            all arrays must be big enough for return
  //            rtn[0] is ntri, rtn[1] is nbound, rtn[2] is ncoord
  //            rtn[3] is nvert, rtn[4] is nfree, rtn[5] is nuniqueb
  //            ntri triangles, t1,t2,t3 are index into xc,yc
  //            nbound boundary edges, b1,b2 index into xc,yc
  //            ncoord vertices or nodes, index order, start at zero
  tri_split(int ntri, int t1[], int t2[], int t3[],
            int nbound, int b1[], int b2[],
            int ncoord, double xc[], double yc[],
            int uniqueb[], int rtn[])
  {
    int j12, j23, j31;
    int nvert, nntri, nnbound, nfree, nuniqueb;
    int i1, i2, i3;
    int nboundt;
    double x1, y1, x2, y2, x3, y3;

    nvert=ncoord;
    nntri=ntri;
    nnbound=nbound;
    nboundt=nbound;

    for(int i=0; i<ntri; i++)
    {
      i1=t1[i];
      i2=t2[i];
      i3=t3[i];
      x1=xc[i1];
      y1=yc[i1];
      x2=xc[i2];
      y2=yc[i2];
      x3=xc[i3];
      y3=yc[i3];

      // create 3 new coordinates, side mid points
      j12=nvert;
      nvert++;
      xc[j12]=(x1+x2)/2.0;
      yc[j12]=(y1+y2)/2.0;
      for(int k=0; k<nvert-1; k++) // do not duplicate
      {
	if(Math.abs(xc[j12]-xc[k])+Math.abs(yc[j12]-yc[k])<0.0001)
	{
	  j12 = k;
          nvert--;
          break;
	}
      }
      j23=nvert;
      nvert++;
      xc[j23]=(x2+x3)/2.0;
      yc[j23]=(y2+y3)/2.0;
      for(int k=0; k<nvert-1; k++) // do not duplicate
      {
	if(Math.abs(xc[j23]-xc[k])+Math.abs(yc[j23]-yc[k])<0.0001)
	{
	  j23 = k;
          nvert--;
          break;
	}
      }
      j31=nvert;
      nvert++;
      xc[j31]=(x3+x1)/2.0;
      yc[j31]=(y3+y1)/2.0;
      for(int k=0; k<nvert-1; k++) // do not duplicate
      {
	if(Math.abs(xc[j31]-xc[k])+Math.abs(yc[j31]-yc[k])<0.0001)
	{
	  j31 = k;
          nvert--;
          break;
	}
      }

      // create 3 new triangles, center replaces original
      t1[i]=j12; t2[i]=j23; t3[i]=j31;
      t1[nntri]=i1; t2[nntri]=j12; t3[nntri]=j31;
      nntri++;
      t1[nntri]=i2; t2[nntri]=j23; t3[nntri]=j12;
      nntri++;
      t1[nntri]=i3; t2[nntri]=j31; t3[nntri]=j23;
      nntri++;

      // if any side is a boundary, add and fix
      for(int j=0; j<nboundt; j++)
      {
	if((b1[j]==i1 && b2[j]==i2) || (b1[j]==i2 && b2[j]==i1))
	{
	  b2[nnbound]=b2[j];
          b1[nnbound]=j12;
          b2[j]=j12;
          nnbound++;
	  break;
	}
      } // end j<nbound
      nboundt = nnbound;
      for(int j=0; j<nboundt; j++)
      {
        if((b1[j]==i1 && b2[j]==i3) || (b1[j]==i3 && b2[j]==i1))
	{
	  b2[nnbound]=b2[j];
          b1[nnbound]=j31;
          b2[j]=j31;
          nnbound++;
	  break;
	}
      } // end j<nbound
      nboundt = nnbound;
      for(int j=0; j<nboundt; j++)
      {
        if((b1[j]==i2 && b2[j]==i3) || (b1[j]==i3 && b2[j]==i2))
	{
	  b2[nnbound]=b2[j];
          b1[nnbound]=j23;
          b2[j]=j23;
          nnbound++;
	  break;
	}
      } // end j<nbound
      nboundt = nnbound;

    } // end i<ntri

    // compute  uniqueb  unique list of boundary nodes
    nuniqueb = 0;
    for(int i=0; i<nnbound; i++)
    {
      uniqueb[nuniqueb]=b1[i];
      uniqueb[nuniqueb+1]=b2[i];
      nuniqueb=nuniqueb+2;
    }
    nuniqueb=uniqueint(nuniqueb, uniqueb);
    nfree = nvert - nuniqueb;
    // caller must store new values
    rtn[0] = nntri;     // ntri
    rtn[1] = nnbound;   // nbound
    rtn[2] = nvert;     // ncoord
    rtn[3] = nvert;     // nvert
    rtn[4] = nfree;     // nfree
    rtn[5] = nuniqueb;  // nuniqueb
  } // end tri_split 

  int uniqueint(int nold, int a[])
  {
    int i, j;
    sortint(nold,a);
    j=0;
    for(i=1; i<nold; i++)
    {
      if(a[i]>a[i-1])
      {
	j++;
	a[j]=a[i];
      }
    }
    j++;
    return j;
  } // end uniqueint

  void sortint(int n, int a[])
  {
    int i, j, temp;
    for(i=0; i<n-1; i++)
      for(j=i+1; j<n; j++)
	if(a[i]>a[j]) {temp=a[i]; a[i]=a[j]; a[j]=temp;}
  } // end sortint
} // end class tri_split