#include #include #include #include #include #include using namespace std; #include #include #include #include #include #include #include #include #include typedef point4 pt4; typedef point4 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 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 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(keval(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=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 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(k13 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) ); }