Back to index

salome-med  6.5.0
SplitterTetra.hxx
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 #ifndef __SPLITTERTETRA_HXX__
00021 #define __SPLITTERTETRA_HXX__
00022 
00023 #include "TransformedTriangle.hxx"
00024 #include "TetraAffineTransform.hxx"
00025 #include "InterpolationOptions.hxx"
00026 #include "InterpKernelHashMap.hxx"
00027 #include "VectorUtils.hxx"
00028 
00029 #include <assert.h>
00030 #include <vector>
00031 #include <functional>
00032 #include <map>
00033 #include <set>
00034 
00035 namespace INTERP_KERNEL
00036 {
00041   class TriangleFaceKey
00042   {
00043   public:
00044     
00054     TriangleFaceKey(int node1, int node2, int node3)
00055     {
00056       sort3Ints(_nodes, node1, node2, node3);
00057       _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
00058     }
00059 
00067     bool operator==(const TriangleFaceKey& key) const
00068     {
00069       return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
00070     }
00071 
00078     bool operator<(const TriangleFaceKey& key) const
00079     {
00080       for (int i = 0; i < 3; ++i)
00081         {
00082           if (_nodes[i] < key._nodes[i])
00083             {
00084               return true;
00085             }
00086           else if (_nodes[i] > key._nodes[i])
00087             {
00088               return false;
00089             }
00090         }
00091       return false;
00092     }
00093 
00100     int hashVal() const
00101     {
00102       return _hashVal;
00103     }
00104      
00105     inline void sort3Ints(int* sorted, int node1, int node2, int node3);
00106 
00107   private:
00109     int _nodes[3];
00110     
00112     int _hashVal;
00113   };
00114   
00123   inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
00124   {
00125     if(x1 < x2)
00126       {
00127         if(x1 < x3)
00128           {
00129             // x1 is min
00130             sorted[0] = x1;
00131             sorted[1] = x2 < x3 ? x2 : x3;
00132             sorted[2] = x2 < x3 ? x3 : x2;
00133           }
00134         else
00135           {
00136             // x3, x1, x2
00137             sorted[0] = x3;
00138             sorted[1] = x1;
00139             sorted[2] = x2;
00140           }
00141       }
00142     else // x2 < x1
00143       {
00144         if(x2 < x3)
00145           {
00146             // x2 is min
00147             sorted[0] = x2;
00148             sorted[1] = x1 < x3 ? x1 : x3;
00149             sorted[2] = x1 < x3 ? x3 : x1;
00150           } 
00151         else
00152           {
00153             // x3, x2, x1
00154             sorted[0] = x3;
00155             sorted[1] = x2;
00156             sorted[2] = x1;
00157           }
00158       }
00159   }
00160 
00166   template<> class hash<INTERP_KERNEL::TriangleFaceKey>
00167   {
00168   public:
00175     int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
00176     {
00177       return key.hashVal();
00178     }
00179   };
00180 }
00181 
00182 namespace INTERP_KERNEL
00183 {
00184 
00190   template<class MyMeshType>
00191   class SplitterTetra
00192   {
00193   public: 
00194     
00195     SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);
00196 
00197     ~SplitterTetra();
00198 
00199     double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
00200     double intersectSourceFace(const NormalizedCellType polyType,
00201                                const int polyNodesNbr,
00202                                const int *const polyNodes,
00203                                const double *const *const polyCoords,
00204                                const double dimCaracteristic,
00205                                const double precision,
00206                                std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
00207                                std::set<TriangleFaceKey>& listOfTetraFacesColinear);
00208 
00209     double intersectTetra(const double** tetraCorners);
00210 
00211     typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
00212     
00213     void splitIntoDualCells(SplitterTetra<MyMeshType> **output);
00214 
00215     void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);
00216 
00217     void clearVolumesCache();
00218 
00219   private:
00220     // member functions
00221     inline void createAffineTransform(const double** corners);
00222     inline void checkIsOutside(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const;
00223     inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
00224     inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
00225     inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
00226     inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
00227     inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
00228 
00229     static inline bool IsFacesCoplanar(const double *const planeNormal, const double planeConstant,
00230                                 const double *const *const coordsFace, const double precision);
00231     static inline double CalculateIntersectionSurfaceOfCoplanarTriangles(const double *const planeNormal,
00232                                                                          const double planeConstant,
00233                                                                          const double *const p1, const double *const p2, const double *const p3,
00234                                                                          const double *const p4, const double *const p5, const double *const p6,
00235                                                                          const double dimCaracteristic, const double precision);
00236 
00238     SplitterTetra(const SplitterTetra& t);
00239     
00241     SplitterTetra& operator=(const SplitterTetra& t);
00242 
00243     // member variables
00245     TetraAffineTransform* _t;
00246     
00248     HashMap< int , double* > _nodes;
00249     
00251     HashMap< TriangleFaceKey, double
00252 // #ifdef WIN32
00253 //         , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator> 
00254 // #endif
00255     > _volumes;
00256 
00258     const MyMeshType& _src_mesh;
00259                 
00260     // node id of the first node in target mesh in C mode.
00261     typename MyMeshType::MyConnType _conn[4];
00262 
00263     double _coords[12];
00264   };
00265 
00272   template<class MyMeshType>
00273   inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
00274   {
00275     // create AffineTransform from tetrahedron
00276     _t = new TetraAffineTransform( corners );
00277   }
00278 
00288   template<class MyMeshType>
00289   inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside, const double errTol) const
00290   {
00291     isOutside[0] = isOutside[0] && (pt[0] < errTol);
00292     isOutside[1] = isOutside[1] && (pt[0] > (1.0-errTol) );
00293     isOutside[2] = isOutside[2] && (pt[1] < errTol);
00294     isOutside[3] = isOutside[3] && (pt[1] > (1.0-errTol));
00295     isOutside[4] = isOutside[4] && (pt[2] < errTol);
00296     isOutside[5] = isOutside[5] && (pt[2] > (1.0-errTol));
00297     isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < errTol);
00298     isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0-errTol) );
00299   }
00300   
00301   template<class MyMeshType>
00302   inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
00303   {
00304     isStrictlyOutside[0] = isStrictlyOutside[0] && (pt[0] < -errTol);
00305     isStrictlyOutside[1] = isStrictlyOutside[1] && (pt[0] > (1.0 + errTol));
00306     isStrictlyOutside[2] = isStrictlyOutside[2] && (pt[1] < -errTol);
00307     isStrictlyOutside[3] = isStrictlyOutside[3] && (pt[1] > (1.0 + errTol));
00308     isStrictlyOutside[4] = isStrictlyOutside[4] && (pt[2] < -errTol);
00309     isStrictlyOutside[5] = isStrictlyOutside[5] && (pt[2] > (1.0 + errTol));
00310     isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
00311     isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
00312   }
00313 
00323   template<class MyMeshType>
00324   inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
00325   {  
00326     const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
00327     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
00328     assert(transformedNode != 0);
00329     _t->apply(transformedNode, node);
00330     _nodes[globalNodeNum] = transformedNode;
00331   }
00332 
00333 
00344   template<class MyMeshType>
00345   inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
00346   {
00347     double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
00348     assert(transformedNode != 0);
00349     _t->apply(transformedNode, node);
00350     _nodes[globalNodeNum] = transformedNode;
00351   }
00352 
00360   template<class MyMeshType>
00361   inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
00362   {
00363     const double vol = tri.calculateIntersectionVolume();
00364     _volumes.insert(std::make_pair(key, vol));
00365   }
00366 
00374   template<class MyMeshType>
00375   inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
00376   {
00377     const double surf = tri.calculateIntersectionSurface(_t);
00378     _volumes.insert(std::make_pair(key, surf));
00379   }
00380 
00381   template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
00382   class SplitterTetra2
00383   {
00384   public:
00385     SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
00386     ~SplitterTetra2();
00387     void releaseArrays();
00388     void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
00389                          typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00390     void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00391     void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00392     void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00393     void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00394     void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00395     void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
00396                      typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
00397     void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
00398     inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
00399     inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
00400     //template<int n>
00401     inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
00402   private:
00403     const MyMeshTypeT& _target_mesh;
00404     const MyMeshTypeS& _src_mesh;
00405     SplittingPolicy _splitting_pol;
00408     std::vector<const double*> _nodes;
00409     std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
00410   };
00411 
00419   template<class MyMeshTypeT, class MyMeshTypeS>
00420   //template<int n>
00421   inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
00422   {
00423     barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
00424     for(int i = 0; i < n ; ++i)
00425       {
00426        const double* pt = getCoordsOfSubNode(pts[i]);
00427        barycenter[0] += pt[0];
00428        barycenter[1] += pt[1];
00429        barycenter[2] += pt[2];
00430       }
00431     
00432     barycenter[0] /= n;
00433     barycenter[1] /= n;
00434     barycenter[2] /= n;
00435   }
00436 
00443   template<class MyMeshTypeT, class MyMeshTypeS>
00444   inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
00445   {
00446     // replace "at()" with [] for unsafe but faster access
00447     return _nodes.at(node);
00448   }
00449 
00457   template<class MyMeshTypeT, class MyMeshTypeS>
00458   const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
00459   {
00460     const double *ret=_nodes.at(node);
00461     if(node<8)
00462       nodeId=_node_ids[node];
00463     else
00464       nodeId=-1;
00465     return ret;    
00466   }
00467 }
00468 
00469 #endif