preCalculateTriangleSurroundsEdge();
preCalculateTripleProducts();
-#if 1 //CS_ALM
_indiceA = 0;
_indiceB = 0;
-#endif
+ kA=0;
+ kB=0;
}
{
// delete elements of polygons
-#if 0 //CS_ALM
- for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
- {
- delete[] *it;
- }
-#else
- //_polygonA.clear();
-#endif
-
-
-
-#if 0 //CS_ALM
- for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
- {
- delete[] *it;
- }
-#else
- //_polygonB.clear();
-
_indiceA = 0;
_indiceB = 0;
-#endif
+ kA=0;
+ kB=0;
}
/**
// calculate volume under A
double volA = 0.0;
-#if 0 //CS_ALM
- if(_polygonA.size() > 2)
-#else
if(_indiceA > 2)
-#endif
{
LOG(2, "---- Treating polygon A ... ");
calculatePolygonBarycenter(A, barycenter);
double volB = 0.0;
// if triangle is not in h = 0 plane, calculate volume under B
-#if 0 //CS_ALM
- if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
-#else
if(_indiceB > 2 && !isTriangleInPlaneOfFacet(XYZ))
-#endif
{
LOG(2, "---- Treating polygon B ... ");
*/
void TransformedTriangle::calculateIntersectionPolygons()
{
-#if 1 //CS_ALM
assert(_indiceA == 0);
assert(_indiceB == 0);
-#endif
// avoid reallocations in push_back() by pre-allocating enough memory
// we should never have more than 20 points
-#if 0 //CS_ALM
- _polygonA.reserve(20);
- _polygonB.reserve(20);
-#else
- //_polygonA.reserve(60);
- //_polygonB.reserve(60);
-
-#endif
// -- surface intersections
// surface - edge
-
-
for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
{
if(testSurfaceEdgeIntersection(edge))
{
-#if 0 //CS_ALM
- double* ptA = new double[3];
- calcIntersectionPtSurfaceEdge(edge, ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
-#else
double ptA[3] = {0.0,0.0,0.0};
calcIntersectionPtSurfaceEdge(edge, ptA);
+
for (int i = 0; i < 3; i++) {
- _polygonA[_indiceA][i] = ptA[i];
+ _polygonA[kA++] = ptA[i];
}
+
_indiceA += 1;
- LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
-#endif
+ LOG(1,"Surface-edge : " << vToStr(ptA) << " added to A ");
if(edge >= XY)
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(ptA, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
-#else
- double ptB[3]= {0.0,0.0,0.0};
- copyVector3(ptA, ptB);
for (int i = 0; i < 3; i++) {
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = ptA[i];
}
+
_indiceB += 1;
- LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
-#endif
+ LOG(1,"Surface-edge : " << vToStr(ptA) << " added to B ");
}
{
if(testSurfaceRayIntersection(corner))
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
-#else
- double ptB[3]= {0.0,0.0,0.0} ;
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
for (int i = 0; i < 3; i++) {
+ _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i];
- _polygonB[_indiceB][i] = ptB[i];
}
+
_indiceB += 1;
- LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
-#endif
+ //LOG(1,"Surface-ray : " << vToStr(pt) << " added to B");
}
if(doTest && testSegmentFacetIntersection(seg, facet))
{
-#if 0 //CS_ALM
- double* ptA = new double[3];
- calcIntersectionPtSegmentFacet(seg, facet, ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
-#else
double ptA[3] = {0.0,0.0,0.0};
calcIntersectionPtSegmentFacet(seg, facet, ptA);
for (int i = 0; i < 3; i++) {
- _polygonA[_indiceA][i] = ptA[i];
+ _polygonA[kA++] = ptA[i];
}
+
_indiceA += 1;
- LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
-#endif
+ LOG(1,"Segment-facet : " << vToStr(ptA) << " added to A");
if(facet == XYZ)
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(ptA, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
-#else
- double ptB[3]= {0.0,0.0,0.0};
- copyVector3(ptA, ptB);
for (int i = 0; i < 3; i++) {
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = ptA[i];
}
+
_indiceB += 1;
- LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
-#endif
+ LOG(1,"Segment-facet : " << vToStr(ptA) << " added to B");
}
if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
{
-#if 0 //CS_ALM
- double* ptA = new double[3];
- calcIntersectionPtSegmentEdge(seg, edge, ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
- if(edge >= XY)
- {
- double* ptB = new double[3];
- copyVector3(ptA, ptB);
- _polygonB.push_back(ptB);
- }
-#else
double ptA[3]= {0.0,0.0,0.0} ;
calcIntersectionPtSegmentEdge(seg, edge, ptA);
for (int i = 0; i < 3; i++) {
- _polygonA[_indiceA][i] = ptA[i];
+ _polygonA[kA++] = ptA[i];
}
+
_indiceA += 1;
- LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
+ LOG(1,"Segment-edge : " << vToStr(ptA) << " added to A");
if(edge >= XY)
{
- double ptB[3]= {0.0,0.0,0.0};
- copyVector3(ptA, ptB);
for (int i = 0; i < 3; i++) {
-
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = ptA[i];
}
+
_indiceB += 1;
}
-#endif
}
}
if(doTest && testSegmentCornerIntersection(seg, corner))
{
-#if 0 //CS_ALM
- double* ptA = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
- if(corner != O)
- {
- double* ptB = new double[3];
- _polygonB.push_back(ptB);
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
- }
-#else
- double ptA[3]= {0.0,0.0,0.0};
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
for (int i = 0; i < 3; i++) {
+ _polygonA[kA++] = COORDS_TET_CORNER[3 * corner+i];
- _polygonA[_indiceA][i] = ptA[i];
}
_indiceA += 1;
- LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
+
+ LOG(1,"Segment-corner : " << vToStr(ptA) << " added to A");
if(corner != O)
{
- double ptB[3]= {0.0,0.0,0.0};
-
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
for (int i = 0; i < 3; i++) {
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i];
+
}
+
_indiceB += 1;
- LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
+ LOG(1,"Segment-corner : " << vToStr(ptB) << " added to B");
}
-#endif
}
}
{
if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
-#else
- double ptB[3]= {0.0,0.0,0.0};
- copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
for (int i = 0; i < 3; i++) {
-
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = COORDS_TET_CORNER[3 * corner+i];
+
}
+
_indiceB += 1;
- LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
-#endif
+ LOG(1,"Segment-ray : " << vToStr(ptB) << " added to B");
}
}
#endif
if(testSegmentHalfstripIntersection(seg, edge))
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
-#else
double ptB[3]= {0.0,0.0,0.0};
calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
for (int i = 0; i < 3; i++) {
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = ptB[i];
}
+
_indiceB += 1;
- LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
-#endif
+ LOG(1,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
}
}
}
// tetrahedron
if(testCornerInTetrahedron(corner))
{
-#if 0 //CS_ALM
- double* ptA = new double[3];
- copyVector3(&_coords[5*corner], ptA);
- _polygonA.push_back(ptA);
- LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
-#else
- double ptA[3]= {0.0,0.0,0.0};
- copyVector3(&_coords[5*corner], ptA);
- for (int i = 0; i < 3; i++) {
- _polygonA[_indiceA][i] = ptA[i];
+ for (int i = 0; i < 3; i++) {
+ _polygonA[kA++] = _coords[5*corner+i];
+
}
+
_indiceA += 1;
- LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
-#endif
- }
+ LOG(1,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
+ }
// XYZ - plane
if(testCornerOnXYZFacet(corner))
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- _polygonB.push_back(ptB);
- LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
-#else
- double ptB[3]= {0.0,0.0,0.0};
- copyVector3(&_coords[5*corner], ptB);
- for (int i = 0; i < 3; i++) {
-
- _polygonB[_indiceB][i] = ptB[i];
+ for (int i = 0; i < 3; i++) {
+ _polygonB[kB++] = _coords[5*corner+i];
+
}
+
_indiceB += 1;
- LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
-#endif
+ LOG(1,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
}
// projection on XYZ - facet
if(testCornerAboveXYZFacet(corner))
{
-#if 0 //CS_ALM
- double* ptB = new double[3];
- copyVector3(&_coords[5*corner], ptB);
- ptB[2] = 1 - ptB[0] - ptB[1];
- assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
- _polygonB.push_back(ptB);
- LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
-#else
double ptB[3]= {0.0,0.0,0.0};
copyVector3(&_coords[5*corner], ptB);
ptB[2] = 1 - ptB[0] - ptB[1];
assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
for (int i = 0; i < 3; i++) {
- _polygonB[_indiceB][i] = ptB[i];
+ _polygonB[kB++] = ptB[i];
}
+
_indiceB += 1;
- LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
-#endif
+ LOG(1,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
}
}
+ for (int i = 0; i < _indiceA; i++) {
+ _polygonA2[i] = &_polygonA[i*3];
+ }
+
+ for (int i = 0; i < _indiceB; i++) {
+ _polygonB2[i] = &_polygonB[i*3];
+ }
+
}
/**
* @param barycenter array of three doubles where barycenter is stored
*
*/
-#if 0 //CS_ALM
void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
-#else
- void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
-#endif
{
LOG(3,"--- Calculating polygon barycenter");
// get the polygon points
-#if 0 //CS_ALM
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
-#else
- VEC3 polygon[20];
+ double polygon[20*3];
+ double *polygon2[20];
if (poly == A) {
- for (int i=0; i < _indiceA; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonA[i][j];
- }
+ for (int i=0; i < _indiceA*3; i++) {
+ polygon[i] = _polygonA[i];
+
}
+ for (int i=0; i < _indiceA; i++) {
+ polygon2[i] = &polygon[i*3];
+ }
}
else {
- for (int i=0; i < _indiceB; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonB[i][j];
- }
+ for (int i=0; i < _indiceB*3; i++) {
+ polygon[i] = _polygonB[i];
+
}
- }
-
-#endif
+ for (int i=0; i < _indiceB; i++) {
+ polygon2[i] = &polygon[i*3];
+ }
+ }
+
+
// calculate barycenter
-#if 0 //CS_ALM
- const int m = polygon.size();
-#else
const int m = (poly == A) ? _indiceA : _indiceB ;
-#endif
-
for(int j = 0 ; j < 3 ; ++j)
{
barycenter[j] = 0.0;
{
for(int i = 0 ; i < m ; ++i)
{
-
- const double* pt = polygon[i];
- for(int j = 0 ; j < 3 ; ++j)
- {
- barycenter[j] += pt[j] / double(m);
- }
+
+ const double* pt = polygon2[i];
+ for(int j = 0 ; j < 3 ; ++j)
+ {
+ barycenter[j] += pt[j] / double(m);
+ }
+
}
}
LOG(3,"Barycenter is " << vToStr(barycenter));
typedef SortOrder::CoordType CoordType;
// get the polygon points
-#if 0 //CS_ALM
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
- if(polygon.size() == 0)
- return;
-#else
- VEC3 polygon[20];
+ double* polygon2[20];
+
if (poly == A) {
+ /*for (int i=0; i < _indiceA*3; i++) {
+ polygon[i] = _polygonA[i];
+
+ }*/
for (int i=0; i < _indiceA; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonA[i][j];
- }
- }
+ polygon2[i] = new double[3];
+ polygon2[i][0] = _polygonA[i*3];
+ polygon2[i][1] = _polygonA[i*3 + 1];
+ polygon2[i][2] = _polygonA[i*3 + 2];
+ // polygon2[i] = &_polygonA[i*3];
+ }
}
else {
+ /*for (int i=0; i < _indiceB*3; i++) {
+
+ polygon[i] = _polygonB[i];
+ }*/
for (int i=0; i < _indiceB; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonB[i][j];
- }
+ // polygon2[i] = &_polygonB[i*3];
+ polygon2[i] = new double[3];
+ polygon2[i][0] = _polygonB[i*3];
+ polygon2[i][1] = _polygonB[i*3 + 1];
+ polygon2[i][2] = _polygonB[i*3 + 2];
}
}
int taillePoly = (poly == A) ? _indiceA : _indiceB;
if (taillePoly == 0)
return;
-#endif
-
// determine type of sorting
CoordType type = SortOrder::XY;
{
type = SortOrder::YZ;
}
- }
+ }
// create order object
SortOrder order(barycenter, type);
// sort vector with this object
// NB : do not change place of first object, with respect to which the order
// is defined
-#if 0 //CS_ALM
- sort((polygon.begin()), polygon.end(), order);
-#else
- //double* pol;
- //std::sort( pol, pol+60, order);
- std::sort( polygon, polygon+(taillePoly*3), order);
-#endif
+ std::sort( polygon2, polygon2+taillePoly, order);
+ if (poly == A) {
+ for (int i=0; i < _indiceA; i++) {
+ // _polygonA[i*3] = *polygon2[i];
+ _polygonA[i*3] = polygon2[i][0];
+ _polygonA[i*3 + 1] = polygon2[i][1];
+ _polygonA[i*3 + 2] = polygon2[i][2];
+
+ }
+ }
+ else {
+ for (int i=0; i < _indiceB; i++) {
+ // _polygonB[i*3] = *polygon2[i];
+ _polygonB[i*3] = polygon2[i][0];
+ _polygonB[i*3 + 1] = polygon2[i][1];
+ _polygonB[i*3 + 2] = polygon2[i][2];
+
+ }
+ }
+
LOG(3,"Sorted polygon is ");
-#if 0 //CS_ALM
- for(size_t i = 0 ; i < polygon.size() ; ++i)
-#else
for(size_t i = 0 ; i < (size_t)taillePoly ; ++i)
-#endif
+
{
- LOG(3,vToStr(polygon[i]));
+ if (poly == A) {
+ LOG(1,vToStr(&_polygonA[i*3]));
+ }
+ else {
+ LOG(1,vToStr(&_polygonB[i*3]));
+ }
}
+// for (int i = 0; i < taillePoly; i++) {
+// delete [] polygon2[i];
+// }
+// delete polygon2;
}
/**
LOG(2,"--- Calculating volume under polygon");
// get the polygon points
-#if 0 //CS_ALM
- std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
- int taillePoly = polygon.size()
-#else
- VEC3 polygon[20];
- if (poly == A) {
+
+ double *polygon2[20];
+ if (poly == A) {
+
for (int i=0; i < _indiceA; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonA[i][j];
- }
- }
+ polygon2[i] = &_polygonA[i*3];
+ }
}
else {
+
for (int i=0; i < _indiceB; i++) {
- for (int j = 0; j < 3; j++) {
- polygon[i][j] = _polygonB[i][j];
- }
- }
- }
+ polygon2[i] = &_polygonB[i*3];
+ }
+ }
int taillePoly = (poly == A) ? _indiceA : _indiceB;
-
-#endif
double vol = 0.0;
const int m = taillePoly;
for(int i = 0 ; i < m ; ++i)
{
-
- const double* ptCurr = polygon[i]; // pt "i"
- const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
-
+ const double* ptCurr = polygon2[i]; // pt "i"
+ const double* ptNext = polygon2[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
const double factor2 =
ptCurr[0]*(ptNext[1] - barycenter[1])
+ barycenter[0]*(ptCurr[1] - ptNext[1]);
vol += (factor1 * factor2) / 6.0;
}
+
LOG(2,"Abs. Volume is " << vol);
return vol;
double triNormal[3], polyNormal[3];
crossprod<3>( triangle.getCorner(P),triangle.getCorner(Q),triangle.getCorner(R), triNormal);
-#if 0 //CS_ALM
- const std::vector<double*> * pPolygonA = &triangle.getPolygonA();
-#else
- const std::vector<double> * pPolygonA = &triangle.getPolygonA();
-#endif
- if ( pPolygonA->size() < 3 )
+
+ const double* pPolygonA = triangle.getPolygonA();
+
+
+ if ( triangle.getIndiceA() < 3 )
{
if ( !epsilonEqual( triNormal[_Z], 0 ))
return; // not vertical triangle does not intersect the unit tetra
// Vertical triangle. Use inherited methods of TransformedTriangle to
// calculate intesection polygon
*((TransformedTriangle*)this) = triangle; // copy triangle fields
- _polygonA.clear();
- _polygonB.clear();
+ this->_indiceA = 0;
+ this->_indiceB= 0;
+ this->kA = 0 ;
+ this->kB = 0;
+
calculateIntersectionPolygons();
- if (this->_polygonA.size() < 3)
+ if (this->_indiceA < 3)
return;
calculatePolygonBarycenter(A, _barycenterA);
sortIntersectionPolygon(A, _barycenterA);
- pPolygonA = & _polygonA;
+
+ pPolygonA = _polygonA;
+
}
+
// check if polygon orientation is same as the one of triangle
-#if 0 //CS_ALM
- std::vector<double*>::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end();
-#else
- std::vector<double>::const_iterator p = pPolygonA->begin(), pEnd = pPolygonA->end();
- /*int taillePoly = pPolygonA->size();
- double p[3];
- double pEnd[3];
-
- for (int j = 0; j < 3; j++){
- p[j] = pPolygonA->operator[](j);
- }*/
-
-#endif
- //#ifdef DMP_UNITTETRAINTERSECTIONBARY
+ for (int i = 0; i < _indiceA; i ++)
+ _polygonA2[i] = &_polygonA[i*3];
+
+ for (int i = 0; i < _indiceB; i ++)
+ _polygonB2[i] = &_polygonB[i*3];
+
+ const double* p = _polygonA2[0];
+ const double* pEnd = _polygonA2[_indiceA];
+
+
+#ifdef DMP_UNITTETRAINTERSECTIONBARY
std::cout.precision(18);
std::cout << "**** int polygon() " << std::endl;
while ( p != pEnd )
p += 3;
}
p = pPolygonA->begin();
- //#endif
+#endif
-#if 0 //CS_ALM
- double* p1 = *p++;
- double* p2 = *p;
- while ( samePoint( p1, p2 ) && ++p != pEnd )
- p2 = *p;
-#else
-
-#endif
+ const double* p1 = p++;
+ const double* p2 = p;
+ while ( samePoint( p1, p2 ) && ++p != pEnd )
+ p2 = p;
if ( p == pEnd )
{
clearPolygons();
return;
}
- double* p3 = *p;
+
+ const double* p3 = p;
while (( samePoint( p2, p3 ) || samePoint( p1, p3 )) && ++p != pEnd )
- p3 = *p;
+ p3 = p;
+
if ( p == pEnd )
{
#ifdef DMP_UNITTETRAINTERSECTIONBARY
_faces.push_back( std::vector< double* > () );
std::vector< double* >& faceCorner = _faces.back();
- faceCorner.resize( pPolygonA->size()/* + 1*/ );
+ faceCorner.resize( _indiceA/* + 1*/ );
int i = 0;
if ( reverse )
{
- std::vector<double*>::const_reverse_iterator polyF = pPolygonA->rbegin(), polyEnd;
- for ( polyEnd = pPolygonA->rend(); polyF != polyEnd; ++i, ++polyF )
- if ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) {
-#if 0 //CS_ALM
- copyVector3( *polyF, faceCorner[i] = new double[3] );
-#else
+
+ const double* polyF = _polygonA2[_indiceA];
+ const double * polyEnd;
+ //for ( polyEnd = _polygonA2[0]; polyF != polyEnd; ++i, ++polyF )
+ for ( polyEnd = _polygonA2[0]; polyF != polyEnd; ++i, --polyF )
+ if ( i==0 || !samePoint( polyF, faceCorner[i-1] )) {
+
+
+
double var[3]= {0.0,0.0,0.0};
- faceCorner[i] = var;
- copyVector3( *polyF, faceCorner[i] );
- //faceCorner[i] = var;
-#endif
+ copyVector3(faceCorner[i] ,var);
+ copyVector3( polyF, faceCorner[i] );
+
+
} else {
--i;
}
}
else
{
- std::vector<double*>::const_iterator polyF = pPolygonA->begin(), polyEnd;
- for ( polyEnd = pPolygonA->end(); polyF != polyEnd; ++i, ++polyF )
- if ( i==0 || !samePoint( *polyF, faceCorner[i-1] )) {
-#if 0 //CS_ALM
- copyVector3( *polyF, faceCorner[i] = new double[3] );
-#else
+
+ //const double* polyF = pPolygonA;
+ const double* polyF = _polygonA2[0];
+ const double * polyEnd;
+ for ( polyEnd = _polygonA2[_indiceA]; polyF != polyEnd; ++i, ++polyF )
+ if ( i==0 || !samePoint( polyF, faceCorner[i-1] )) {
+
double var[3]= {0.0,0.0,0.0};
- faceCorner[i] = var;
- copyVector3( *polyF, faceCorner[i] );
-#endif
+ copyVector3(faceCorner[i] ,var);
+ copyVector3( polyF, faceCorner[i] );
+
} else {
--i;
}
if ( i < 3 )
{
clearPolygons(); // free memory of _polygonA
- _polygonA = faceCorner;
+
+ int facesize = faceCorner.size();
+ int k = 0;
+ for (int i = 0; i < facesize*3; i++,k++) {
+ for (int j = 0; j < 3; j++) {
+ _polygonA[k++] = faceCorner[i][j];
+ }
+
+ }
+
_faces.pop_back();
}
else
{
- if ( i < (int)pPolygonA->size() )
+
+ if ( i < (int)_indiceA )
faceCorner.resize( i );
if ( _polyNormals.empty() )
double bary[] = { 0, 0, 0 };
// base polygon bary
-#if 1 //CS_ALM
+
int polySize = polygon.size();
-#endif
+
for ( v = polygon.begin(); v != vEnd ; ++v )
{
double* p = *v;
bary[1] += p[1];
bary[2] += p[2];
}
-#if 0 //CS_ALM
- bary[0] /= polygon.size();
- bary[1] /= polygon.size();
- bary[2] /= polygon.size();
-#else
+
bary[0] /= polySize;
bary[1] /= polySize;
bary[2] /= polySize;
-#endif
+
// pyramid volume
double vol = 0;
-#if 0 //CS_ALM
- for ( int i = 0; i < (int)polygon.size(); ++i )
-#else
+
for ( int i = 0; i < polySize; ++i )
-#endif
+
{
double* p1 = polygon[i];
-#if 0 //CS_ALM
- double* p2 = polygon[(i+1)%polygon.size()];
-#else
+
double* p2 = polygon[(i+1)%polySize];
-#endif
+
vol += std::fabs( calculateVolumeForTetra( p1, p2, bary, P ));
}
{
std::vector< double* >& polygon = *f;
double coordSum[3] = {0,0,0};
-#if 1 //CS_ALM
+
int polySize = polygon.size();
-#endif
-#if 0 //CS_ALM
- for ( int i = 0; i < (int)polygon.size(); ++i )
-#else
for ( int i = 0; i < polySize; ++i )
-#endif
{
double* p = polygon[i];
coordSum[0] += p[0];
if ( epsilonEqual( coordSum[j], 0.0 ))
sideAdded[j] = ++nbAddedSides != 0 ;
}
-#if 0 //CS_ALM
- if ( !sideAdded[3] &&
- ( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polygon.size(), 1. ))) {
-#else
+
if ( !sideAdded[3] &&
( epsilonEqual( (coordSum[0]+coordSum[1]+coordSum[2]) / polySize, 1. ))) {
-#endif
sideAdded[3] = ++nbAddedSides != 0 ;
}
}
// ---------------------------------------------------------------------------------
int nbIntersectPolygs = _faces.size();
-#if 1 //CS_ALM
+
std::vector< double* > * sideFaces[ 4 ]; // future polygons on sides of tetra
-#else
- std::vector< double[3] > * sideFaces[ 4 ]; // future polygons on sides of tetra
-#endif
+
for ( int i = 0; i < NB_TETRA_SIDES; ++i )
{
sideFaces[ i ]=0;
for ( int iF = 0; iF < nbIntersectPolygs; ++f, ++iF ) // loop on added intersection polygons
{
std::vector< double* >& polygon = *f;
-#if 1 //CS_ALM
+
int polySize = (int)polygon.size();
-#endif
+
for ( int i = 0; i < polySize; ++i )
{
// segment ends
double* p1 = polygon[i];
-#if 0 //CS_ALM
- double* p2 = polygon[(i+1)%polygon.size()];
-#else
+
double* p2 = polygon[(i+1)%polySize];
-#endif
+
bool p1OnSide, p2OnSide;//, onZeroSide = false;
for ( int j = 0; j < 3; ++j )
{
if ( p1OnSide && p2OnSide )
{
// segment p1-p2 is on j-th orthogonal side of tetra
-#if 0 //CS_ALM
- sideFaces[j]->push_back( new double[3] );
- copyVector3( p1, sideFaces[j]->back() );
- sideFaces[j]->push_back( new double[3] );
- copyVector3( p2, sideFaces[j]->back() );
-#else
+
double var[3]= {0.0,0.0,0.0};
sideFaces[j]->push_back( var );
copyVector3( p1, sideFaces[j]->back() );
double var2[3]= {0.0,0.0,0.0};
sideFaces[j]->push_back( var2 );
copyVector3( p2, sideFaces[j]->back() );
-#endif
+
//break;
}
epsilonEqual( p1[_X] + p1[_Y] + p1[_Z], 1.0 ) &&
epsilonEqual( p2[_X] + p2[_Y] + p2[_Z], 1.0 ))
{
-#if 0 //CS_ALM
- sideFaces[3]->push_back( new double[3] );
- copyVector3( p1, sideFaces[3]->back() );
- sideFaces[3]->push_back( new double[3] );
- copyVector3( p2, sideFaces[3]->back() );
-#else
+
double var[3]= {0.0,0.0,0.0};
sideFaces[3]->push_back( var );
copyVector3( p1, sideFaces[3]->back() );
double var2[3]= {0.0,0.0,0.0};
sideFaces[3]->push_back( var2 );
copyVector3( p2, sideFaces[3]->back() );
-#endif
+
}
}
}
else
{
std::vector< double* >& sideFace = *sideFaces[i];
-#if 0 //CS_ALM
- for ( int i = 0; i < sideFace.size(); ++i )
- {
-#else
+
int faceSize = sideFace.size();
for ( int i = 0; i < faceSize; ++i )
{
-#endif
double* p = sideFace[i];
std::cout << "\t" << i << ": ( " << p[0] << ", " << p[1] << ", " << p[2] << " )" << std::endl;
{
if ( !cutOffCorners[ ic ] && ic != excludeCorner )
{
-#if 0 //CS_ALM
- sideFace.push_back( new double[3] );
-#else
+
double var[3]= {0.0,0.0,0.0};
sideFace.push_back( var );
-#endif
+
copyVector3( origin, sideFace.back() );
if ( ic )
sideFace.back()[ ic-1 ] = 1.0;
if ( face.size() >= 3 )
{
clearPolygons(); // free memory of _polygonA
- _polygonA = face;
+
+ int k = 0;
+ for (unsigned int i = 0; i < face.size(); i++,k++) {
+ for (int j = 0; j < 3; j++) {
+ _polygonA[k++] = face[i][j];
+ }
+ _indiceA++;
+ }
face.clear();
- face.reserve( _polygonA.size() );
+ face.reserve( _indiceA );
+
+
if ( iF >= nbIntersectPolygs )
{ // sort points of side faces
calculatePolygonBarycenter( A, _barycenterA );
sortIntersectionPolygon( A, _barycenterA );
}
// exclude equal points
- std::vector< double* >::iterator v = _polygonA.begin(), vEnd = _polygonA.end();
- face.push_back( *v );
- *v = 0;
- for ( ++v; v != vEnd; ++v )
+
+ for (int i = 0; i < _indiceA; i++) {
+ _polygonA2[i] = &_polygonA[i*3];
+ }
+ double* v = _polygonA2[0];
+ double *vEnd = _polygonA2[_indiceA];
+ face.push_back( v );
+
+
+
+
+ *v = 0;
+ for ( ++v; v != vEnd; ++v )
{
- double* pPrev = face.back();
- double* p = *v;
+ double* pPrev = face.back();
+ double* p = v;
if ( !samePoint( p, pPrev ))
- {
- face.push_back( p );
- *v = 0;
- }
- }
+ {
+ face.push_back( p );
+ *v = 0;
+ }
+ }
+
}
if ( face.size() < 3 )
{ // size could decrease
clearPolygons(); // free memory of _polygonA
- _polygonA = face;
- face.clear();
+
+ int k = 0;
+ for (unsigned int i = 0; i < face.size(); i++,k++) {
+ for (int j = 0; j < 3; j++) {
+ _polygonA[k++] = face[i][j];
+ }
+ }
+ face.clear();
+
}
else
{
void UnitTetraIntersectionBary::clearPolygons(bool andFaces)
{
-#if 0 //CS_ALM
- for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
- {
- delete[] *it;
- *it = 0;
- }
-
- for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
- {
-
- delete[] *it;
-
- *it = 0;
- }
-#endif
-
- _polygonA.clear();
- _polygonB.clear();
+ _indiceA=0;
+ _indiceB=0;
+ kA=0;
+ kB=0;
- if ( andFaces )
- {
- std::list< std::vector< double* > >::iterator f = this->_faces.begin(), fEnd = this->_faces.end();
-#if 0 //CS_ALM
- for ( ; f != fEnd; ++f )
- {
- std::vector< double* >& polygon = *f;
- for(std::vector<double*>::iterator it = polygon.begin() ; it != polygon.end() ; ++it)
- {
-
- delete[] *it;
- *it = 0;
- }
- }
-#endif
- this->_faces.clear();
- }
}
//================================================================================
UnitTetraIntersectionBary::~UnitTetraIntersectionBary()
{
- clearPolygons(/*andFaces=*/true );
+ //clearPolygons(/*andFaces=*/true );
}
}
