Back to index

salome-med  6.5.0
TransformedTriangle.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "TransformedTriangle.hxx"
00021 #include "VectorUtils.hxx"
00022 #include "TetraAffineTransform.hxx"
00023 #include <iostream>
00024 #include <fstream>
00025 #include <cassert>
00026 #include <algorithm>
00027 #include <functional>
00028 #include <iterator>
00029 #include <math.h>
00030 #include <vector>
00031 
00032 namespace INTERP_KERNEL
00033 {
00034 
00040   class ProjectedCentralCircularSortOrder
00041   {
00042   public:
00043 
00045     enum CoordType { XY, XZ, YZ };
00046   
00054     ProjectedCentralCircularSortOrder(const double* barycenter, const CoordType type)
00055       : _aIdx((type == YZ) ? 1 : 0), 
00056         _bIdx((type == XY) ? 1 : 2),
00057         _a(barycenter[_aIdx]), 
00058         _b(barycenter[_bIdx])
00059     {
00060     }
00061 
00071     bool operator()(const double* pt1, const double* pt2)
00072     {
00073       // calculate angles with the axis
00074       const double ang1 = atan2(pt1[_aIdx] - _a, pt1[_bIdx] - _b);
00075       const double ang2 = atan2(pt2[_aIdx] - _a, pt2[_bIdx] - _b);
00076 
00077       return ang1 > ang2;
00078     }
00079 
00080   private:
00082     const int _aIdx;
00083   
00085     const int _bIdx;
00086 
00088     const double _a;
00089   
00091     const double _b;
00092   };
00093 
00094   // ----------------------------------------------------------------------------------
00095   // TransformedTriangle PUBLIC  
00096   // ----------------------------------------------------------------------------------
00097   
00107   TransformedTriangle::TransformedTriangle(double* p, double* q, double* r)
00108     : _is_double_products_calculated(false),  _is_triple_products_calculated(false), _volume(0)
00109   {
00110   
00111     for(int i = 0 ; i < 3 ; ++i) 
00112       {
00113         // xyz coordinates
00114         _coords[5*P + i] = p[i];
00115         _coords[5*Q + i] = q[i];
00116         _coords[5*R + i] = r[i];
00117       }
00118 
00119     // h coordinate
00120     
00121     _coords[5*P + 3] = 1 - p[0] - p[1] - p[2];
00122     _coords[5*Q + 3] = 1 - q[0] - q[1] - q[2];
00123     _coords[5*R + 3] = 1 - r[0] - r[1] - r[2];
00124 
00125     // H coordinate
00126     _coords[5*P + 4] = 1 - p[0] - p[1];
00127     _coords[5*Q + 4] = 1 - q[0] - q[1];
00128     _coords[5*R + 4] = 1 - r[0] - r[1];
00129 
00130     resetNearZeroCoordinates();
00131 
00132     // initialise rest of data
00133     preCalculateDoubleProducts();
00134 
00135     preCalculateTriangleSurroundsEdge();
00136 
00137     preCalculateTripleProducts();
00138  
00139   }
00140 
00147   TransformedTriangle::~TransformedTriangle()
00148   {
00149     // delete elements of polygons
00150     for(std::vector<double*>::iterator it = _polygonA.begin() ; it != _polygonA.end() ; ++it)
00151       {
00152         delete[] *it;
00153       }
00154     for(std::vector<double*>::iterator it = _polygonB.begin() ; it != _polygonB.end() ; ++it)
00155       {
00156         delete[] *it;
00157       }    
00158   }
00159 
00168   double TransformedTriangle::calculateIntersectionVolume()
00169   {
00170     // check first that we are not below z - plane    
00171     if(isTriangleBelowTetraeder())
00172       {
00173         LOG(2, " --- Triangle is below tetraeder - V = 0.0");
00174         return 0.0;
00175       }
00176     
00177     // get the sign of the volume -  equal to the sign of the z-component of the normal
00178     // of the triangle, u_x * v_y - u_y * v_x, where u = q - p and v = r - p
00179     // if it is zero, the triangle is perpendicular to the z - plane and so the volume is zero
00180 //     const double uv_xy[4] = 
00181 //       {
00182 //         _coords[5*Q] - _coords[5*P], _coords[5*Q + 1] - _coords[5*P + 1], // u_x, u_y
00183 //         _coords[5*R] - _coords[5*P], _coords[5*R + 1] - _coords[5*P + 1]  // v_x, v_y
00184 //       };
00185 
00186 //     double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
00187     int sign = isTriangleInclinedToFacet( OXY );
00188 
00189     if(sign == 0 )
00190       {
00191         LOG(2, " --- Triangle is perpendicular to z-plane - V = 0.0");
00192         return _volume = 0.0;
00193       }
00194 
00195 
00196     // normalize sign
00197     //sign = sign > 0.0 ? 1.0 : -1.0;
00198 
00199     LOG(2, "-- Calculating intersection polygons ... ");
00200     calculateIntersectionAndProjectionPolygons();
00201     
00202     double barycenter[3];
00203 
00204     // calculate volume under A
00205     double volA = 0.0;
00206     if(_polygonA.size() > 2)
00207       {
00208         LOG(2, "---- Treating polygon A ... ");
00209         calculatePolygonBarycenter(A, barycenter);
00210         sortIntersectionPolygon(A, barycenter);
00211         volA = calculateVolumeUnderPolygon(A, barycenter);
00212         LOG(2, "Volume is " << sign * volA);
00213       }
00214 
00215     double volB = 0.0;
00216     // if triangle is not in h = 0 plane, calculate volume under B
00217     if(_polygonB.size() > 2 && !isTriangleInPlaneOfFacet(XYZ))
00218       {
00219         LOG(2, "---- Treating polygon B ... ");
00220        
00221         calculatePolygonBarycenter(B, barycenter);
00222         sortIntersectionPolygon(B, barycenter);
00223         volB = calculateVolumeUnderPolygon(B, barycenter);
00224         LOG(2, "Volume is " << sign * volB);
00225       }
00226 
00227     LOG(2, "volA + volB = " << sign * (volA + volB) << std::endl << "***********");
00228   
00229     return _volume = sign * (volA + volB);
00230 
00231   } 
00232 
00241   double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat)
00242   {
00243     // check first that we are not below z - plane
00244     if(isTriangleBelowTetraeder())
00245       {
00246         LOG(2, " --- Triangle is below tetraeder - V = 0.0");
00247         return 0.0;
00248       }
00249 
00250     LOG(2, "-- Calculating intersection polygon ... ");
00251     calculateIntersectionPolygon();
00252 
00253     _volume = 0.;
00254     if(_polygonA.size() > 2) {
00255       double barycenter[3];
00256       calculatePolygonBarycenter(A, barycenter);
00257       sortIntersectionPolygon(A, barycenter);
00258       const std::size_t nbPoints = _polygonA.size();
00259       for(std::size_t i = 0 ; i < nbPoints ; ++i)
00260         tat->reverseApply(_polygonA[i], _polygonA[i]);
00261       _volume = calculateSurfacePolygon();
00262     }
00263 
00264     return _volume;
00265   }
00266 
00267   // ----------------------------------------------------------------------------------
00268   // TransformedTriangle PRIVATE
00269   // ----------------------------------------------------------------------------------
00270 
00279   void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
00280   {
00281     assert(_polygonA.size() == 0);
00282     assert(_polygonB.size() == 0);
00283     // avoid reallocations in push_back() by pre-allocating enough memory
00284     // we should never have more than 20 points
00285     _polygonA.reserve(20);
00286     _polygonB.reserve(20);
00287     // -- surface intersections
00288     // surface - edge
00289     for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
00290       {
00291         if(testSurfaceEdgeIntersection(edge))
00292           {
00293             double* ptA = new double[3];
00294             calcIntersectionPtSurfaceEdge(edge, ptA);
00295             _polygonA.push_back(ptA);
00296             LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
00297             if(edge >= XY)
00298               {
00299                 double* ptB = new double[3];
00300                 copyVector3(ptA, ptB);
00301                 _polygonB.push_back(ptB);
00302                 LOG(3,"Surface-edge : " << vToStr(ptB) << " added to B ");
00303               }
00304            
00305           }
00306       }
00307 
00308     // surface - ray
00309     for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
00310       {
00311         if(testSurfaceRayIntersection(corner))
00312           {
00313             double* ptB = new double[3];
00314             copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
00315             _polygonB.push_back(ptB);
00316             LOG(3,"Surface-ray : " << vToStr(ptB) << " added to B");
00317           }
00318       }
00319 
00320     // -- segment intersections
00321     for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
00322       {
00323 
00324         bool isZero[NO_DP];
00325 
00326         // check beforehand which double-products are zero
00327         for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
00328           {
00329             isZero[dp] = (calcStableC(seg, dp) == 0.0);
00330           }
00331 
00332         // segment - facet
00333         for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
00334           {
00335             // is this test worth it?
00336             const bool doTest = 
00337               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] && 
00338               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
00339               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
00340 
00341             if(doTest && testSegmentFacetIntersection(seg, facet))
00342               {
00343                 double* ptA = new double[3];
00344                 calcIntersectionPtSegmentFacet(seg, facet, ptA);
00345                 _polygonA.push_back(ptA);
00346                 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
00347                 if(facet == XYZ)
00348                   {
00349                     double* ptB = new double[3];
00350                     copyVector3(ptA, ptB);
00351                     _polygonB.push_back(ptB);
00352                     LOG(3,"Segment-facet : " << vToStr(ptB) << " added to B");
00353                   }
00354 
00355               }
00356           }
00357 
00358         // segment - edge
00359         for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
00360           {
00361             const DoubleProduct edge_dp = DoubleProduct(edge);
00362 
00363             if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
00364               {
00365                 double* ptA = new double[3];
00366                 calcIntersectionPtSegmentEdge(seg, edge, ptA);
00367                 _polygonA.push_back(ptA);
00368                 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
00369                 if(edge >= XY)
00370                   {
00371                     double* ptB = new double[3];
00372                     copyVector3(ptA, ptB);
00373                     _polygonB.push_back(ptB);
00374                   }
00375               }
00376           }
00377        
00378         // segment - corner
00379         for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
00380           {
00381             const bool doTest = 
00382               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
00383               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
00384               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
00385 
00386             if(doTest && testSegmentCornerIntersection(seg, corner))
00387               {
00388                 double* ptA = new double[3];
00389                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
00390                 _polygonA.push_back(ptA);
00391                 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
00392                 if(corner != O)
00393                   {
00394                     double* ptB = new double[3];
00395                     _polygonB.push_back(ptB);
00396                     copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
00397                     LOG(3,"Segment-corner : " << vToStr(ptB) << " added to B");
00398                   }
00399               }
00400           }
00401 
00402             // segment - ray 
00403             for(TetraCorner corner = X ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
00404               {
00405                 if(isZero[DP_SEGMENT_RAY_INTERSECTION[7*(corner-1)]] && testSegmentRayIntersection(seg, corner))
00406                   {
00407                     double* ptB = new double[3];
00408                     copyVector3(&COORDS_TET_CORNER[3 * corner], ptB);
00409                     _polygonB.push_back(ptB);
00410                     LOG(3,"Segment-ray : " << vToStr(ptB) << " added to B");
00411                   }
00412               }
00413        
00414             // segment - halfstrip
00415             for(TetraEdge edge = XY ; edge <= ZX ; edge = TetraEdge(edge + 1))
00416               {
00417 
00418 #if 0
00419                 const int edgeIdx = int(edge) - 3; // offset since we only care for edges XY - ZX
00420                 const bool doTest = 
00421                   !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx]] &&
00422                   !isZero[DP_FOR_HALFSTRIP_INTERSECTION[4*edgeIdx+1]];
00423        
00424 
00425                 if(doTest && testSegmentHalfstripIntersection(seg, edge))
00426 #endif
00427                   if(testSegmentHalfstripIntersection(seg, edge))
00428                     {
00429                       double* ptB = new double[3];
00430                       calcIntersectionPtSegmentHalfstrip(seg, edge, ptB);
00431                       _polygonB.push_back(ptB);
00432                       LOG(3,"Segment-halfstrip : " << vToStr(ptB) << " added to B");
00433                     }
00434               }
00435       }
00436 
00437         // inclusion tests
00438         for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
00439           {
00440             // { XYZ - inclusion only possible if in Tetrahedron?
00441             // tetrahedron
00442             if(testCornerInTetrahedron(corner))
00443               {
00444                 double* ptA = new double[3];
00445                 copyVector3(&_coords[5*corner], ptA);
00446                 _polygonA.push_back(ptA);
00447                 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
00448               }
00449 
00450             // XYZ - plane
00451             if(testCornerOnXYZFacet(corner))
00452               {
00453                 double* ptB = new double[3];
00454                 copyVector3(&_coords[5*corner], ptB);
00455                 _polygonB.push_back(ptB);
00456                 LOG(3,"Inclusion XYZ-plane : " << vToStr(ptB) << " added to B");
00457               }
00458 
00459             // projection on XYZ - facet
00460             if(testCornerAboveXYZFacet(corner))
00461               {
00462                 double* ptB = new double[3];
00463                 copyVector3(&_coords[5*corner], ptB);
00464                 ptB[2] = 1 - ptB[0] - ptB[1];
00465                 assert(epsilonEqual(ptB[0]+ptB[1]+ptB[2] - 1, 0.0));
00466                 _polygonB.push_back(ptB);
00467                 LOG(3,"Projection XYZ-plane : " << vToStr(ptB) << " added to B");
00468               }
00469 
00470           }
00471 
00472       }
00473 
00481   void TransformedTriangle::calculateIntersectionPolygon()
00482   {
00483     assert(_polygonA.size() == 0);
00484     // avoid reallocations in push_back() by pre-allocating enough memory
00485     // we should never have more than 20 points
00486     _polygonA.reserve(20);
00487     // -- surface intersections
00488     // surface - edge
00489     for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
00490       {
00491         if(testSurfaceEdgeIntersection(edge))
00492           {
00493             double* ptA = new double[3];
00494             calcIntersectionPtSurfaceEdge(edge, ptA);
00495             _polygonA.push_back(ptA);
00496             LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A ");
00497           }
00498       }
00499 
00500     // -- segment intersections
00501     for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1))
00502       {
00503 
00504         bool isZero[NO_DP];
00505 
00506         // check beforehand which double-products are zero
00507         for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1))
00508           {
00509             isZero[dp] = (calcStableC(seg, dp) == 0.0);
00510           }
00511 
00512         // segment - facet
00513         for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1))
00514           {
00515             // is this test worth it?
00516             const bool doTest =
00517               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] &&
00518               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] &&
00519               !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]];
00520 
00521             if(doTest && testSegmentFacetIntersection(seg, facet))
00522               {
00523                 double* ptA = new double[3];
00524                 calcIntersectionPtSegmentFacet(seg, facet, ptA);
00525                 _polygonA.push_back(ptA);
00526                 LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A");
00527               }
00528           }
00529 
00530         // segment - edge
00531         for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1))
00532           {
00533             const DoubleProduct edge_dp = DoubleProduct(edge);
00534 
00535             if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge))
00536               {
00537                 double* ptA = new double[3];
00538                 calcIntersectionPtSegmentEdge(seg, edge, ptA);
00539                 _polygonA.push_back(ptA);
00540                 LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A");
00541               }
00542           }
00543 
00544         // segment - corner
00545         for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1))
00546           {
00547             const bool doTest =
00548               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] &&
00549               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] &&
00550               isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )];
00551 
00552             if(doTest && testSegmentCornerIntersection(seg, corner))
00553               {
00554                 double* ptA = new double[3];
00555                 copyVector3(&COORDS_TET_CORNER[3 * corner], ptA);
00556                 _polygonA.push_back(ptA);
00557                 LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A");
00558               }
00559           }
00560 
00561       }
00562 
00563         // inclusion tests
00564         for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1))
00565           {
00566             // { XYZ - inclusion only possible if in Tetrahedron?
00567             // tetrahedron
00568             if(testCornerInTetrahedron(corner))
00569               {
00570                 double* ptA = new double[3];
00571                 copyVector3(&_coords[5*corner], ptA);
00572                 _polygonA.push_back(ptA);
00573                 LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A");
00574               }
00575 
00576           }
00577 
00578       }
00579 
00580 
00586     double TransformedTriangle::calculateSurfacePolygon()
00587     {
00588       const std::size_t nbPoints = _polygonA.size();
00589       double pdt[3];
00590       double sum[3] = {0., 0., 0.};
00591 
00592       for(std::size_t i = 0 ; i < nbPoints ; ++i)
00593         {
00594           const double *const ptCurr = _polygonA[i];  // pt "i"
00595           const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt nbPoints == pt 0)
00596 
00597           cross(ptCurr, ptNext, pdt);
00598           add(pdt, sum);
00599         }
00600 
00601       const double surface = norm(sum) * 0.5;
00602       LOG(2,"Surface is " << surface);
00603       return surface;
00604     }
00605 
00615     void TransformedTriangle::calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter)
00616     {
00617       LOG(3,"--- Calculating polygon barycenter");
00618 
00619       // get the polygon points
00620       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
00621 
00622       // calculate barycenter
00623       const std::size_t m = polygon.size();
00624 
00625       for(int j = 0 ; j < 3 ; ++j)
00626         {
00627           barycenter[j] = 0.0;
00628         }
00629 
00630       if(m != 0)
00631         {
00632           for(std::size_t i = 0 ; i < m ; ++i)
00633             {
00634               const double* pt = polygon[i];
00635               for(int j = 0 ; j < 3 ; ++j)
00636                 {
00637                   barycenter[j] += pt[j] / double(m);
00638                 }
00639             }
00640         }
00641       LOG(3,"Barycenter is " << vToStr(barycenter));
00642     }
00643 
00654     void TransformedTriangle::sortIntersectionPolygon(const IntersectionPolygon poly, const double* barycenter)
00655     {
00656       LOG(3,"--- Sorting polygon ...");
00657 
00658       using INTERP_KERNEL::ProjectedCentralCircularSortOrder;
00659       typedef ProjectedCentralCircularSortOrder SortOrder; // change is only necessary here and in constructor
00660       typedef SortOrder::CoordType CoordType;
00661 
00662       // get the polygon points
00663       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
00664 
00665       if(polygon.size() == 0)
00666         return;
00667 
00668       // determine type of sorting
00669       CoordType type = SortOrder::XY;
00670       if(poly == A && !isTriangleInclinedToFacet( OXY )) // B is on h = 0 plane -> ok
00671         {
00672           // NB : the following test is never true if we have eliminated the
00673           // triangles parallel to x == 0 and y == 0 in calculateIntersectionVolume().
00674           // We keep the test here anyway, to avoid interdependency.
00675 
00676           // is triangle inclined to x == 0 ?
00677           if(isTriangleInclinedToFacet(OZX))
00678             {
00679               type = SortOrder::XZ;
00680             }
00681           else //if(isTriangleParallelToFacet(OYZ))
00682             {
00683               type = SortOrder::YZ;
00684             }
00685         }
00686 
00687       // create order object
00688       SortOrder order(barycenter, type);
00689 
00690       // sort vector with this object
00691       // NB : do not change place of first object, with respect to which the order
00692       // is defined
00693       sort((polygon.begin()), polygon.end(), order);
00694     
00695       LOG(3,"Sorted polygon is ");
00696       for(size_t i = 0 ; i < polygon.size() ; ++i)
00697         {
00698           LOG(3,vToStr(polygon[i]));
00699         }
00700 
00701     }
00702 
00714     double TransformedTriangle::calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter)
00715     {
00716       LOG(2,"--- Calculating volume under polygon");
00717 
00718       // get the polygon points
00719       std::vector<double*>& polygon = (poly == A) ? _polygonA : _polygonB;
00720 
00721       double vol = 0.0;
00722       const std::size_t m = polygon.size();
00723 
00724       for(std::size_t i = 0 ; i < m ; ++i)
00725         {
00726           const double* ptCurr = polygon[i];  // pt "i"
00727           const double* ptNext = polygon[(i + 1) % m]; // pt "i+1" (pt m == pt 0)
00728        
00729           const double factor1 = ptCurr[2] + ptNext[2] + barycenter[2];
00730           const double factor2 = 
00731             ptCurr[0]*(ptNext[1] - barycenter[1]) 
00732             + ptNext[0]*(barycenter[1] - ptCurr[1])
00733             + barycenter[0]*(ptCurr[1] - ptNext[1]);
00734           vol += (factor1 * factor2) / 6.0;
00735         }
00736 
00737       LOG(2,"Abs. Volume is " << vol); 
00738       return vol;
00739     }
00740 
00741 
00743     // Detection of (very) degenerate cases                                /////////////
00745 
00752     bool TransformedTriangle::isTriangleInPlaneOfFacet(const TetraFacet facet) const
00753     {
00754 
00755       // coordinate to check
00756       const int coord = static_cast<int>(facet);
00757 
00758       for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
00759         {
00760           if(_coords[5*c + coord] != 0.0)
00761             {
00762               return false;
00763             }
00764         }
00765     
00766       return true;
00767     }
00768 
00775     bool TransformedTriangle::isTriangleParallelToFacet(const TetraFacet facet) const
00776     {
00777       // coordinate to check
00778       const int coord = static_cast<int>(facet);
00779       return (_coords[5*P + coord] == _coords[5*Q + coord]) && (_coords[5*P + coord] == _coords[5*R + coord]);
00780     }
00781 
00789     int TransformedTriangle::isTriangleInclinedToFacet(const TetraFacet facet) const
00790     {
00791       // coordinate to check
00792       const int coord = static_cast<int>(facet);
00793       const int ind1 = ( coord+1 ) % 3, ind2 = ( coord+2 ) % 3;
00794       const double uv_xy[4] = 
00795         {
00796           // u_x, u_y
00797           _coords[5*Q+ind1] - _coords[5*P+ind1], _coords[5*Q+ind2] - _coords[5*P+ind2],
00798           // v_x, v_y
00799           _coords[5*R+ind1] - _coords[5*P+ind1], _coords[5*R+ind2] - _coords[5*P+ind2]
00800         };
00801 
00802       double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
00803       if(epsilonEqual(sign, 0.))
00804         {
00805           sign = 0.;
00806         }
00807       return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
00808     }
00809 
00815     bool TransformedTriangle::isTriangleBelowTetraeder() const
00816     {
00817       for(TriCorner c = P ; c < NO_TRI_CORNER ; c = TriCorner(c + 1))
00818         {
00819           // check z-coords for all points
00820           if(_coords[5*c + 2] >= 0.0)
00821             {
00822               return false;
00823             }
00824         }
00825       return true;
00826     }
00827 
00832     void TransformedTriangle::dumpCoords() const
00833     {
00834       std::cout << "Coords : ";
00835       for(int i = 0 ; i < 3; ++i)
00836         {
00837           std::cout << vToStr(&_coords[5*i]) << ",";
00838         }
00839       std::cout << std::endl;
00840     }
00841     
00842   } // NAMESPACE