d4tess.cpp

Go to the documentation of this file.

#include <cassert>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
using namespace std;

#include <point.h>
#include <print.h>

#include <d4tri.h>
#include <d4tess.h>

#include <tetrahedron.h>
#include <tetrahedronpartition.h>
#include <halfspaceD3.h>

#include <d4fan.h>
#include <message.h>

typedef point4<double> pt4;
typedef point4<double> const pt4c;


void d4tess::getpoints
( 
  pt4 & P0, 
  pt4 & P1, 
  pt4 & P2, 
  pt4 & P3 
) const
{
  d4tri const & x(cpsimplex());

//messageglobal() << "cp=" << cp << "\n";
//messageglobal() << "x=" << x << "\n";

  P0 = pt[ x.pi[ vs.v[0] ] ];
  P1 = pt[ x.pi[ vs.v[1] ] ];
  P2 = pt[ x.pi[ vs.v[2] ] ];
  P3 = pt[ x.pi[ vs.v[3] ] ];
}


d4tess::~d4tess()
{
  delete minimizer;
  minimizer=0;
}

d4tess::d4tess( uintc arrayMax )
  : fan(*this), minimizer( new d4minoperator(*this) )
{ 
  pt.reserve(arrayMax);
  vi.reserve(arrayMax);

  reset();
}

void d4tess::minimizerSet( d4minoperator* m)
{
  delete minimizer;
  minimizer = m;
}


void d4tess::viadd
(
  uintc v0, uintc v1, uintc v2, uintc v3,
  uintc n0, uintc n1, uintc n2, uintc n3
)
{
  vi.push_back( d4tri(v0,v1,v2,v3,n0,n1,n2,n3) );
}



void d4tess::addpoint(uintc k)
{
  uint face;
  bool found = searchinsidemesh(face,k);

  if (found)
    split(k);
  else
    fan.eval(k);

}











boolc d4tess::surfaceviewable(uintc w) const
{
  //d3halfspace<double> h
  halfspaceD3< pt4, double > h
  (
    pt[ vi[cp].pi[ vs.v[0] ] ],
    pt[ vi[cp].pi[ vs.v[1] ] ],
    pt[ vi[cp].pi[ vs.v[2] ] ]
  );

assert(false);  // changed halfspace code.
  return h.isInside(pt[w]); 
//  return h.eval(pt[w]);
}








void d4tess::infiniteRecursion()
{
  messageglobal() << "\n"; 
  messageglobal() << "*************************************************" << "\n";
  messageglobal() << "       Find the error." << "\n";
  messageglobal() << "\n";
  messageglobal() << *this << "\n";
  messageglobal() << "error:  Infinite Recursion has occured, see previous mesh." << "\n";
}

boolc d4tess::searchinsidemesh( uint & viewableface, uintc p )
{
  bool notfound = false;
  bool isinside = false;
  for ( ; !notfound; )
    notfound = move_terminated(isinside,viewableface,p);

  if (isinside==false)
  {
    if (vs.v[3]!=viewableface)
      vs.left();
    else
      return isinside;

    if (vs.v[3]!=viewableface)
      vs.left();
    else
      return isinside;

    if (vs.v[3]!=viewableface)
      vs.right();
    else
      return isinside;

     
    if (vs.v[3]!=viewableface)
      assert(false);
  }
  
  return isinside;
}

boolc d4tess::move_terminated
(
  bool & insidetetrahedron,
  uint & viewableface,
  uintc p
)
{
  d4tri & x(vi[cp]);

  assert(p<pt.size());
  pt4c & target(pt[p]);

  tetrahedronpartition< pt4, double > tp
  ( 
    pt[ x.pi[0] ],
    pt[ x.pi[1] ],
    pt[ x.pi[2] ],
    pt[ x.pi[3] ]
  );

  // Is the point inside the tetrahedron?
  if (tp.isInside(target))
  {
    insidetetrahedron=true;

    return true;
  }

  for ( uint i=0; i<4; ++i )
  {
    // Is the target viewable?
    if (tp.hi[i].isInside(target))
    {
      if(x.ni[i]==0)
      {
        viewableface=i;
        return true;
      }

      // Move
      cp = x.ni[i];

      return false;
      //i=4; // exit loop
    }
  }

  assert(false);

  return false;
}


void d4tess::debugprinttet(uintc k) const 
{
  uint i;
  for (i=0; i<4; ++i)
    messageglobal() << vi[k].ni[i] << " ";
  messageglobal() << "\n";

  
  for (i=0; i<4; ++i)
    messageglobal() << vi[k].pi[i] << " ";
  messageglobal() << "\n";
}


// Four tetrahedrons from one.
void d4tess::split( uintc k )
{
  assert(k<pt.size());
  assert(k!=0);

  uintc w=cp;


/*
  uintc n0 = vi[w].ni[0];
  uintc n1 = vi[w].ni[1];
  uintc n2 = vi[w].ni[2];
  uintc n3 = vi[w].ni[3];
*/

  uint n[4];
  uint p[4];

  for (uint i=0; i<4; ++i)
  {
    n[i] = vi[w].ni[i];
    p[i] = vi[w].pi[i];
  }

/*

messageglobal() << SHOW(n0) << " ";
messageglobal() << SHOW(n1) << " ";
messageglobal() << SHOW(n2) << " ";
messageglobal() << SHOW(n3) << " "; 
messageglobal() << "\n";

  uintc p0 = vi[w].pi[0];
  uintc p1 = vi[w].pi[1];
  uintc p2 = vi[w].pi[2];
  uintc p3 = vi[w].pi[3];

messageglobal() << SHOW(p0) << " ";
messageglobal() << SHOW(p1) << " ";
messageglobal() << SHOW(p2) << " ";
messageglobal() << SHOW(p3) << " "; 
messageglobal() << "\n";
*/

  // Reference to three new consecutive tetrahedrons.
  // Which are the tetrahedrons opposite 0,1,2 of cp(before) respecitely.
  // The base tetrahedron w will now be the new cp.
  uintc z=vi.size();

  vi[w] =       d4tri(p[0],p[1],p[2],k,  z,z+1,z+2, n[3]);     // w:

  vi.push_back( d4tri(p[3],p[2],p[1],k,  w,z+2,z+1, n[0]) );   // z:
  vi.push_back( d4tri(p[3],p[0],p[2],k,  w,z,z+2,   n[1]) );   // z+1: 
  vi.push_back( d4tri(p[3],p[1],p[0],k,  w,z+1,z,   n[2]) );   // z+2

  if (n[0]!=0)
    vi[n[0]].niReset(cp,z);
  if (n[1]!=0)
    vi[n[1]].niReset(cp,z+1);
  if (n[2]!=0)
    vi[n[2]].niReset(cp,z+2);



    //vi[n0].ni[ vi[n0].niInverse(cp) ] = z;


//  print(messageglobal());

/*
messageglobal() << vi[w] << "\n";
messageglobal() << vi[z] << "\n";
messageglobal() << vi[z+1] << " ";
messageglobal() << vi[z+2] << " ";
*/

/*
  if (n1!=0)
    vi[n1].ni[ vi[n1].niInverse(cp) ] = z+1;
  if (n2!=0)
    vi[n2].ni[ vi[n2].niInverse(cp) ] = z+2;
*/

  bool write_error(false);

  if (write_error)
  {
    messagefile e1("e1.txt",true); 
    e1() << *this;
  }
  debugcheck();

  if (n[0]!=0)
    minimizer->eval(n[0],z);
  if (write_error)
  {
    messagefile e2("e2.txt",true); 
    e2() << *this;
  }
  debugcheck();

  if (n[1]!=0)
    minimizer->eval(n[1],z+1);
  if (write_error)
  {
    messagefile e3("e3.txt",true); 
    e3() << *this;
  }
  debugcheck();

  if (n[2]!=0)
    minimizer->eval(n[2],z+2);
  if (write_error)
  {
    messagefile e4("e4.txt",true); 
    e4() << *this;
  }
  debugcheck();

  if (n[3]!=0)
    minimizer->eval(n[3],w);
  if (write_error)
  {
    messagefile e5("e5.txt",true); 
    e5() << *this;
  }
  debugcheck();
}

ostream & d4tess::print(ostream & os) const
{
  uintc ptsz = pt.size();
  os << SHOW(ptsz) << "\n";

  for (uint i=1; i<ptsz; ++i)
  {
    os << setw(2) << i;
    os << "  " << pt[i] << "\n";
  }

  uintc visz = vi.size();
  os << SHOW(visz) << "\n";
  for (uint i=1; i<visz; ++i)
  {
    os << setw(2) << i;
    os << "  " << vi[i] << "\n";
  }
     
  return os;
}


boolc d4tess::valid(uintc k) const
{
  if (k>=vi.size())
    return false;

  if (vi[k].isnull())
    return true;

  // Look at each neighbour for links.
  for (uint i=0, nb; i<4; ++i)
  {
    nb = vi[k].ni[i];
    if (nb==0)
      continue;

    bool res;

    // Each point bordering neighbour must also be one
    // of the neighbours own points.
    for (uint w=0; w<4; ++w)
    {
      if (w==i)
        continue;

      vi[nb].piInverse( res, vi[k].pi[w] );
      if (res==false)
        return false;
    }

    // The nieghbour must have a pointer back to k.
    vi[nb].niInverse( res, k );
    if (res==false)
      return false;
  }

  return true;
}



boolc d4tess::valid() const
{
  bool res(true);

  uintc visz = vi.size();
  for (uint i=1; i<visz; ++i)
    res &= valid(i);

  return res;
}

boolc d4tess::debugcheck() const 
{

#ifndef NDEBUG

  bool res(true);
  uintc visz = vi.size();

  for (uint i=1; i<visz; ++i )
  {
    if ( vi[i].isnull() )
      continue;
    for (uint k=i+1; k<visz; ++k)
    {
      if (vi[i].hassamepoints( vi[k] ))
      {
        res=false;

        messageglobal() << "\n";
        messageglobal() << "error: tetrahedrons with the same points." << "\n";
        messageglobal() << "  check " << i << " " << k << "\n";
        messageglobal() << "vi[i]=" << vi[i] << "\n";
        messageglobal() << "vi[k]=" << vi[k] << "\n";
        messageglobal() << "\n";
        messageglobal() << *this;
      }
    }
  }
  assert(res==true);

  for (uint i=1; i<visz; ++i)
  {
    if(valid(i)==false)
    {
      res=false;
      messageglobal() << "error:  ";
      //messageglobal() << setw(2) << i;
      messageglobal() << " " << i;
      messageglobal() << "  " << vi[i] << "\n";
      messageglobal() << "Printing out tess." << "\n";
      messageglobal() << *this << "\n";
    }
  }

  assert(res==true);

  return res;

#else
  return true;
#endif

}






boolc d4tess::tetmovedown()
{
  d4tri const & x(vi[cp]);

  uint cp2 = x.ni[ vs.v[2] ];
  if (cp2==0)
    return false;

  d4tri const &  y(vi[cp2]);

  uint a = y.piInverse( x.pi[vs.v[0]] );
  uint b = y.piInverse( x.pi[vs.v[1]] );
  uint c = y.piInverse( x.pi[vs.v[3]] );
  uint d = y.niInverse(cp);

  vs.set(a,b,c,d);

  cp = cp2;

  return true;
}















boolc d4tess::tetmoveright()
{
  d4tri const & x(vi[cp]);
messageglobal() << x << "\n";

messageglobal() << "vs.v[0]" << vs.v[0] << "\n";

  uint cp2 = x.ni[ vs.v[0] ];
messageglobal() << "cp2=" << cp2 << "\n";
  if (cp2==0)
    return false;

  d4tri const &  y(vi[cp2]);

messageglobal() << y << "\n";

  uint a = y.piInverse( x.pi[vs.v[1]] );
  uint b = y.piInverse( x.pi[vs.v[2]] );
  uint c = y.piInverse( x.pi[vs.v[3]] );
  uint d = y.niInverse(cp);


  messageglobal() << "a=" << a << " ";
  messageglobal() << "b=" << b << " ";
  messageglobal() << "c=" << c << " ";
  messageglobal() << "d=" << d << " ";
  messageglobal() << "\n";

  vs.set(a,b,c,d);

  cp = cp2;

  return true;
}

boolc d4tess::tetmoveleft()
{
  d4tri const & x(vi[cp]);
messageglobal() << x << "\n";

  uint cp2 = x.ni[ vs.v[1] ];
messageglobal() << "cp2=" << cp2 << "\n";
  if (cp2==0)
    return false;

  d4tri const &  y(vi[cp2]);

messageglobal() << y << "\n";

  uint a = y.piInverse( x.pi[vs.v[2]] );
  uint b = y.piInverse( x.pi[vs.v[0]] );
  uint c = y.piInverse( x.pi[vs.v[3]] );
  uint d = y.niInverse(cp);


  messageglobal() << "a=" << a << " ";
  messageglobal() << "b=" << b << " ";
  messageglobal() << "c=" << c << " ";
  messageglobal() << "d=" << d << " ";
  messageglobal() << "\n";

  vs.set(a,b,c,d);

  cp = cp2;

  return true;
}

void d4tess::surfacedown()
{
  for ( ; tetmovedown(); );

  vs.right();
  vs.right();
  vs.anticlockwise();
}

void d4tess::surfaceleft()
{
  vs.clockwise();
  surfacedown();
}

void d4tess::surfaceright()
{
  vs.anticlockwise();
  surfacedown();
}

boolc d4tess::boundaryorient()
{
  d4tri const & x(vi[cp]);

  // Transverse the virtual tetrahedrons surface.

  if (x.ni[ vs.v[3] ]==0)
    return true;
    
  vs.left();
  if (x.ni[ vs.v[3] ]==0)
    return true;

  vs.left();
  if (x.ni[ vs.v[3] ]==0)
    return true;

  vs.right();
  if (x.ni[ vs.v[3] ]==0)
    return true;

  return false;
}


ostream & operator << (ostream& os, d4tess const & x)
{ 
  return x.print(os); 
}


boolc d4tess::rotateaboutaxis(uint & steps)
{
  steps=0;

  uint cp0 = cp;

  if (tetmovedown()==false)
    return false; 

  steps=1;

  for(;cp!=cp0; ++steps)
  {
    if (tetmovedown()==false)
      return false;
  }

  return true;
}




boolc d4tess::isconnected(uintc a, uintc b) const 
{
  uintc sz(vi.size());
  assert(a<sz);
  assert(b<sz);

  d4tri const & A(vi[a]);
  d4tri const & B(vi[b]);

  bool res;

  A.niInverse(res,b);
  if (res==false)
    return false;

  B.niInverse(res,a);
  if (res==false)
    return false;

  return true;
}


boolc d4tess::isconvex( uintc a, uintc b ) const
{
  assert(a<vi.size());
  assert(b<vi.size());

  assert(a!=0);
  assert(b!=0);

  d4tri const & A(vi[a]);
  d4tri const & B(vi[b]);

//  if (isconnected(a,b)==false)
//    return false;

//  uint ai = A.niInverse(b);

  bool res;
  uint ai = A.niInverse(res,b);
  if (res==false)
    return false;

  uint bi = B.niInverse(res,a);
  if (res==false)
    return false;

// KEEP THIS CODE   [KEEP]
/*

  //  ASSUMING A CONVEX MESH

  //
  // This test may become more efficent as the
  //   mesh is better balanced.
  // At the moment it captures less than 10% of cases.

  for (uint i=0; i<4; ++i)
  {
    if (i==ai)
      continue;

    B.niInverse(res,A.ni[i]);
    if (res==true)
    {
messageglobal() << "non-convex found with integer arithmetic" << " ";
      return false;
    }
  }

messageglobal() << "normal convex test with halfplanes" << " ";

*/



  tetrahedronpartition< pt4, double > tp
  ( 
    pt[A.pi[0] ],
    pt[A.pi[1] ],
    pt[A.pi[2] ],
    pt[A.pi[3] ]
  );

  pt4 p(pt[B.pi[bi]]);

  for (uint i=0; i<4; ++i)
  {
    if (i==ai)
      continue;

    if ( tp.hi[i].isInside(p)==true )
      return false;
  }

  return true;
}

boolc d4tess::tet2to3()
{
  d4tri & x(vi[cp]);
  uint b = x.ni[vs.v[3]];

  return tet2to3_(cp,b);
}




boolc d4tess::tet2to3_(uintc k1, uintc k2)
{
  uintc z = vi.size();
  assert(k1<z);
  assert(k2<z);

  if ((k1==0)||(k2==0))
    return false;


/*
if (isconvex(k1,k2))
  messageglobal() << "convex " << k1 << " " << k2 << "\n";
*/


  if (isconvex(k1,k2)==false)
    return false;

  uint p[5];
  uint n[6];
  uint u[3];

  d4tri const & K1(vi[k1]);
  d4tri const & K2(vi[k2]);

  uintc k1i = K1.niInverse(k2);
  uintc k2i = K2.niInverse(k1);
  
  p[4] = K2.pi[k2i];
  p[0] = K1.pi[k1i];


  K2.getanticlockwiseface( u[0],u[1],u[2], k2i );
  for (uint i=0; i<3; ++i)
    p[i+1] = K2.pi[u[i]]; 

//  n[3] = K2.ni[ u[2] ];
//  n[4] = K2.ni[ u[0] ];
//  n[5] = K2.ni[ u[1] ];

  n[3] = K2.ni[ K2.piInverse(p[3]) ];
  n[4] = K2.ni[ K2.piInverse(p[1]) ];
  n[5] = K2.ni[ K2.piInverse(p[2]) ];

  n[0] = K1.ni[ K1.piInverse(p[3]) ];
  n[1] = K1.ni[ K1.piInverse(p[1]) ];
  n[2] = K1.ni[ K1.piInverse(p[2]) ];
  
  d4tri a[3];

  a[0].construct( p[4],p[0],p[2],p[1],  n[0],n[3],z+2,z+1 ); 
  a[1].construct( p[0],p[4],p[2],p[3],  n[4],n[1],z+2,z ); 
  //a[2].construct( p[0],p[4],p[1],p[3],  n[5],n[2],z+1,z ); 
  a[2].construct( p[0],p[4],p[3],p[1],  n[5],n[2],z,z+1 ); 

  for (uint i=0; i<3; ++i)
    vi.push_back(a[i]);

  if( n[0]!=0)
    vi[n[0]].niReset(k1,z);
  if( n[3]!=0)
    vi[n[3]].niReset(k2,z);
  if( n[1]!=0)
    vi[n[1]].niReset(k1,z+1);
  if( n[4]!=0)
    vi[n[4]].niReset(k2,z+1);
  if( n[2]!=0)
    vi[n[2]].niReset(k1,z+2);
  if( n[5]!=0)
    vi[n[5]].niReset(k2,z+2);

  vi[k1].setnull();
  vi[k2].setnull();

  cp = z;

  return true;
}


boolc d4tess::tet2to3(uintc a, uintc b)
{
//assert(false);
  assert(a<vi.size());
  assert(b<vi.size());

  //messageglobal() << "a=" << a << " b=" << b << " convex " << isconvex(a,b) << "\n";

  if (isconvex(a,b)==false)
    return false;

  uint aj[5];

  d4tri const A(vi[a]);
  aj[0] = A.niInverse(b);

  d4tri const B(vi[b]);
  aj[4] = B.niInverse(a);

  A.getclockwiseface(aj[1],aj[2],aj[3],aj[0]);

  uint pj[5];
  for (uint i=0; i<4; ++i)
    pj[i] = A.pi[aj[i]];
  pj[4] = B.pi[aj[4]];

/*
for (uint i=0; i<5; ++i)
  messageglobal() << "aj[" << i << "]=" << aj[i] << " ";
messageglobal() << "\n";
for (uint i=0; i<5; ++i)
  messageglobal() << "pj[" << i << "]=" << pj[i] << " ";
messageglobal() << "\n";
*/



  // nj[0] is ignored.
  // the indexes correspond with the points : the neighbour is opposite the point
  // with the same index.  For neighbours >3 simply subtract 3 and that is the
  // index for the second tetrahedron.
  uint nj[7];
  
  for (uint i=1; i<=3; ++i)
    nj[i] = A.ni[aj[i]];

  for (uint i=4; i<=6; ++i)
    nj[i] = B.ni[ B.piInverse(pj[i-3]) ];

//for (uint i=1; i<=6; ++i)
//  messageglobal() << "nj[" << i << "]=" << nj[i] << " ";
//messageglobal() << "\n";

  uint k1(a), k2(b), k3(vi.size());

  vi[a] = d4tri( pj[0],pj[4],pj[2],pj[3], nj[4],nj[1],k3,k2 );
  vi[b] = d4tri( pj[4],pj[0],pj[2],pj[1], nj[3],nj[6],k3,k1 );
  vi.push_back( d4tri( pj[4],pj[0],pj[1],pj[3], nj[2],nj[5],k1,k2 ) );

  if (nj[1]!=0)
    vi[nj[1]].niReset(a,k1);
  if (nj[2]!=0)
    vi[nj[2]].niReset(a,k3);
  if (nj[3]!=0)
    vi[nj[3]].niReset(a,k2);
  if (nj[4]!=0)
    vi[nj[4]].niReset(b,k1);
  if (nj[5]!=0)
    vi[nj[5]].niReset(b,k3);
  if (nj[6]!=0)
    vi[nj[6]].niReset(b,k2);

  debugcheck();

  return true;
}

boolc d4tess::tet2to3Inverse()
{
  uint steps;
  if (rotateaboutaxis(steps)==false)
    return false;

  if (steps!=3)
    return false;

  uint p[5];
  uint a[3];
  uint n[6];

  p[0] = vi[cp].pi[ vs.v[0] ];
  p[4] = vi[cp].pi[ vs.v[1] ];

  for (uint i=0; i<3; ++i)
  {
    a[i] = cp;
    p[i+1] = vi[cp].pi[ vs.v[2] ];

    n[i]   = vi[cp].ni[ vs.v[1] ];
    n[i+3] = vi[cp].ni[ vs.v[0] ];

    tetmovedown();
  }

  messageglobal() << "p[i]: ";
  for (uint i=0; i<5; ++i)
    messageglobal() << p[i] << " ";
  messageglobal() << "\n";

  messageglobal() << "n[i] ";
  for (uint i=0; i<6; ++i)
    messageglobal() << n[i] << " ";
  messageglobal() << "\n";

  messageglobal() << "a[i]" << "\n";
  for (uint i=0; i<3; ++i)
    messageglobal() << a[i] << ": " << vi[a[i]] << "\n";
  messageglobal() << "\n";

  uint k = vi.size();
  d4tri k1(p[1],p[3],p[2],p[0], n[1],n[0],n[2],k+1); 
  d4tri k2(p[1],p[2],p[3],p[4], n[4],n[5],n[3],k); 

  messageglobal() << "k1: " << k1 << "\n";
  messageglobal() << "k2: " << k2 << "\n";

  vi.push_back(k1);
  vi.push_back(k2);

  for (uint i=0; i<3; ++i)
    vi[a[0]] = vi[0];

  for (uint i=0; i<3; ++i)
  {
    if (n[i]!=0)
      vi[n[i]].niReset(a[i],k);

    if (n[i+3]!=0)
      vi[n[i+3]].niReset(a[i],k+1);

    vi[a[i]].setnull();
  }

  cp = k;

  debugcheck();


  // Rotating about an axis with three tetrahedrons
  // is always a convex shape.
  return true;
}

simplexface const d4tess::cpsimplexface() const
{ 
  assert(vs.validstate()); 
  return simplexface(cp,vs.v[3]); 
}

void d4tess::reset()
{
  pt.clear();
  vi.clear();

  pt.push_back( pt4(0.0,0.0,0.0,0.0) );
  vi.push_back( d4tri() );

  cp=1;
}

void d4tess::initialize()
{
  assert(pt.size()>4);
  assert(vi.size()==1);

  pt4c & p0(pt[1]);
  pt4c & p1(pt[2]);
  pt4c & p2(pt[3]);
  pt4c & p3(pt[4]);

  halfspaceD3< pt4, double > h(p2,p1,p0);
  if (h.isInside(p3)==false)
    vi.push_back( d4tri(1,2,3,4, 0,0,0,0) );
  else
    vi.push_back( d4tri(4,3,2,1, 0,0,0,0) );
}