src/INTERP_KERNEL/SplitterTetra.hxx

Go to the documentation of this file.

// Copyright (C) 2007-2011  CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//

#ifndef __SPLITTERTETRA_HXX__
#define __SPLITTERTETRA_HXX__

#include "TransformedTriangle.hxx"
#include "TetraAffineTransform.hxx"
#include "InterpolationOptions.hxx"
#include "InterpKernelHashMap.hxx"

#include <assert.h>
#include <vector>
#include <functional>
#include <map>

namespace INTERP_KERNEL
{
  class TriangleFaceKey
  {
  public:
    
    TriangleFaceKey(int node1, int node2, int node3)
    {
      sort3Ints(_nodes, node1, node2, node3);
      _hashVal = ( _nodes[0] + _nodes[1] + _nodes[2] ) % 29;
    }

    bool operator==(const TriangleFaceKey& key) const
    {
      return _nodes[0] == key._nodes[0] && _nodes[1] == key._nodes[1] && _nodes[2] == key._nodes[2];
    }

    int hashVal() const
    {
      return _hashVal;
    }
     
    inline void sort3Ints(int* sorted, int node1, int node2, int node3);

  private:
    int _nodes[3];
    
    int _hashVal;
  };
  
  inline void TriangleFaceKey::sort3Ints(int* sorted, int x1, int x2, int x3)
  {
    if(x1 < x2)
      {
        if(x1 < x3)
          {
            // x1 is min
            sorted[0] = x1;
            sorted[1] = x2 < x3 ? x2 : x3;
            sorted[2] = x2 < x3 ? x3 : x2;
          }
        else
          {
            // x3, x1, x2
            sorted[0] = x3;
            sorted[1] = x1;
            sorted[2] = x2;
          }
      }
    else // x2 < x1
      {
        if(x2 < x3)
          {
            // x2 is min
            sorted[0] = x2;
            sorted[1] = x1 < x3 ? x1 : x3;
            sorted[2] = x1 < x3 ? x3 : x1;
          } 
        else
          {
            // x3, x2, x1
            sorted[0] = x3;
            sorted[1] = x2;
            sorted[2] = x1;
          }
      }
  }

  template<> class hash<INTERP_KERNEL::TriangleFaceKey>
  {
  public:
    int operator()(const INTERP_KERNEL::TriangleFaceKey& key) const
    {
      return key.hashVal();
    }
  };
}

namespace INTERP_KERNEL
{

  template<class MyMeshType>
  class SplitterTetra
  {
  public: 
    
    SplitterTetra(const MyMeshType& srcMesh, const double** tetraCorners, const typename MyMeshType::MyConnType *nodesId);

    ~SplitterTetra();

    double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);

    double intersectTetra(const double** tetraCorners);

    typename MyMeshType::MyConnType getId(int id) { return _conn[id]; }
    
    void splitIntoDualCells(SplitterTetra<MyMeshType> **output);

    void splitMySelfForDual(double* output, int i, typename MyMeshType::MyConnType& nodeId);

    void clearVolumesCache();

  private:
    // member functions
    inline void createAffineTransform(const double** corners);
    inline void checkIsOutside(const double* pt, bool* isOutside) const;
    inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
    inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
        

    SplitterTetra(const SplitterTetra& t);
    
    SplitterTetra& operator=(const SplitterTetra& t);

    // member variables
    TetraAffineTransform* _t;
    
    HashMap< int , double* > _nodes;
    
    HashMap< TriangleFaceKey, double
// #ifdef WIN32
//         , hash_compare<TriangleFaceKey,TriangleFaceKeyComparator> 
// #endif
    > _volumes;

    const MyMeshType& _src_mesh;
                
    // node id of the first node in target mesh in C mode.
    typename MyMeshType::MyConnType _conn[4];

    double _coords[12];
  };

  template<class MyMeshType>
  inline void SplitterTetra<MyMeshType>::createAffineTransform(const double** corners)
  {
    // create AffineTransform from tetrahedron
    _t = new TetraAffineTransform( corners );
  }

  template<class MyMeshType>
  inline void SplitterTetra<MyMeshType>::checkIsOutside(const double* pt, bool* isOutside) const
  {
    isOutside[0] = isOutside[0] && (pt[0] <= 0.0);
    isOutside[1] = isOutside[1] && (pt[0] >= 1.0);
    isOutside[2] = isOutside[2] && (pt[1] <= 0.0);
    isOutside[3] = isOutside[3] && (pt[1] >= 1.0);
    isOutside[4] = isOutside[4] && (pt[2] <= 0.0);
    isOutside[5] = isOutside[5] && (pt[2] >= 1.0);
    isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] <= 0.0);
    isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
  }
  
  template<class MyMeshType>
  inline void SplitterTetra<MyMeshType>::calculateNode(typename MyMeshType::MyConnType globalNodeNum)
  {  
    const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum;
    double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
    assert(transformedNode != 0);
    _t->apply(transformedNode, node);
    _nodes[globalNodeNum] = transformedNode;
  }

  template<class MyMeshType>
  inline void SplitterTetra<MyMeshType>::calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key)
  {
    const double vol = tri.calculateIntersectionVolume();
    _volumes.insert(std::make_pair(key, vol));
  }

  template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
  class SplitterTetra2
  {
  public:
    SplitterTetra2(const MyMeshTypeT& targetMesh, const MyMeshTypeS& srcMesh, SplittingPolicy policy);
    ~SplitterTetra2();
    void releaseArrays();
    void splitTargetCell(typename MyMeshTypeT::MyConnType targetCell, typename MyMeshTypeT::MyConnType nbOfNodesT,
                         typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void fiveSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void sixSplit(const int* const subZone, typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void calculateGeneral24Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void calculateGeneral48Tetra(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void splitPyram5(typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void splitConvex(typename MyMeshTypeT::MyConnType                     targetCell,
                     typename std::vector< SplitterTetra<MyMeshTypeS>* >& tetra);
    void calculateSubNodes(const MyMeshTypeT& targetMesh, typename MyMeshTypeT::MyConnType targetCell);
    inline const double* getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node);
    inline const double* getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId);
    //template<int n>
    inline void calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts);
  private:
    const MyMeshTypeT& _target_mesh;
    const MyMeshTypeS& _src_mesh;
    SplittingPolicy _splitting_pol;
    std::vector<const double*> _nodes;
    std::vector<typename MyMeshTypeT::MyConnType> _node_ids;
  };

  template<class MyMeshTypeT, class MyMeshTypeS>
  //template<int n>
  inline void SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::calcBarycenter(int n, double* barycenter, const typename MyMeshTypeT::MyConnType* pts)
  {
    barycenter[0] = barycenter[1] = barycenter[2] = 0.0;
    for(int i = 0; i < n ; ++i)
      {
       const double* pt = getCoordsOfSubNode(pts[i]);
       barycenter[0] += pt[0];
       barycenter[1] += pt[1];
       barycenter[2] += pt[2];
      }
    
    barycenter[0] /= n;
    barycenter[1] /= n;
    barycenter[2] /= n;
  }

  template<class MyMeshTypeT, class MyMeshTypeS>
  inline const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode(typename MyMeshTypeT::MyConnType node)
  {
    // replace "at()" with [] for unsafe but faster access
    return _nodes.at(node);
  }

  template<class MyMeshTypeT, class MyMeshTypeS>
  const double* SplitterTetra2<MyMeshTypeT, MyMeshTypeS>::getCoordsOfSubNode2(typename MyMeshTypeT::MyConnType node, typename MyMeshTypeT::MyConnType& nodeId)
  {
    const double *ret=_nodes.at(node);
    if(node<8)
      nodeId=_node_ids[node];
    else
      nodeId=-1;
    return ret;    
  }
}

#endif