src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20 #ifndef __PolyhedronIntersectorP1P0Bary_TXX__
  21 #define __PolyhedronIntersectorP1P0Bary_TXX__
  22
  23 #include "PolyhedronIntersectorP1P0Bary.hxx"
  24 #include "Intersector3DP1P0Bary.txx"
  25 #include "MeshUtils.hxx"
  26
  27 #include "SplitterTetra.txx"
  28
  29 namespace INTERP_KERNEL
  30 {
  31
  32   /**
  33    * Constructor creating object from target cell global number
  34    * The constructor first calculates the necessary nodes,
  35    * (depending on the splitting policy) and then splits the hexahedron into
  36    * tetrahedra, placing these in the internal vector _tetra.
  37    *
  38    * @param targetMesh  mesh containing the target elements
  39    * @param srcMesh     mesh containing the source elements
  40    * @param policy      splitting policy to be used
  41    *
  42    * WARNING : in _split attribute, sourceMesh and targetMesh are switched in order to fit intersectCells feature.
  43    */
  44   template<class MyMeshType, class MyMatrix>
  45   PolyhedronIntersectorP1P0Bary<MyMeshType,MyMatrix>::PolyhedronIntersectorP1P0Bary(const MyMeshType& targetMesh,
  46                                                                                     const MyMeshType& srcMesh,
  47                                                                                     SplittingPolicy policy)
  48     :Intersector3DP1P0Bary<MyMeshType,MyMatrix>(targetMesh,srcMesh),_split(targetMesh,srcMesh,policy)
  49   {
  50     // SPEC:
  51     // "Limitation. For the P1P0 barycentric improvement only triangle source cells in 2D and
  52     // tetrahedrons in 3D will be supported by interpolators. If a non
  53     // triangle/tetrahedron source cell is detected an INTERP_KERNEL::Exception should be thrown."
  54
  55     // Check types of source elements here rather than in intersectCells() since a wrong type can be
  56     // found late after a long time of calculation.
  57
  58     const unsigned long numSrcElems = srcMesh.getNumberOfElements();
  59     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
  60       if ( srcMesh.getTypeOfElement( OTT<ConnType,numPol>::indFC(i) ) != NORM_TETRA4 )
  61         throw INTERP_KERNEL::Exception("P1P0 barycentric algorithm works only with tetrahedral source meshes");
  62   }
  63
  64   /**
  65    * Destructor.
  66    * Liberates the SplitterTetra objects and potential sub-node points that have been allocated.
  67    *
  68    */
  69   template<class MyMeshType, class MyMatrix>
  70   PolyhedronIntersectorP1P0Bary<MyMeshType,MyMatrix>::~PolyhedronIntersectorP1P0Bary()
  71   {
  72     releaseArrays();
  73   }
  74
  75   template<class MyMeshType, class MyMatrix>
  76   void PolyhedronIntersectorP1P0Bary<MyMeshType,MyMatrix>::releaseArrays()
  77   {
  78     for(typename std::vector< SplitterTetra<MyMeshType>* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
  79       delete *iter;
  80     _split.releaseArrays();
  81     _tetra.clear();
  82   }
  83
  84   //================================================================================
  85   /*!
  86    * \brief This method computes a value per each node of source cell for each target cell.
  87    *  \param srcCell - a source tetrahedron
  88    *  \param tgtCells - target elements
  89    *  \param res - matrix to fill in
  90    */
  91   //================================================================================
  92
  93   template<class MyMeshType, class MyMatrix>
  94   void PolyhedronIntersectorP1P0Bary<MyMeshType,MyMatrix>::intersectCells(ConnType                     tgtCell,
  95                                                                           const std::vector<ConnType>& srcCells,
  96                                                                           MyMatrix&                    res)
  97   {
  98     typename MyMatrix::value_type& resRow=res[tgtCell];
  99
 100     int nbOfNodesT=Intersector3D<MyMeshType,MyMatrix>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(tgtCell));
 101     releaseArrays();
 102     _split.splitTargetCell(tgtCell,nbOfNodesT,_tetra);
 103
 104     for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
 105     {
 106       // intersect a source tetrahedron with each target tetrahedron: get intersection volume and barycenter
 107       double baryCentre[SPACEDIM], total_baryCentre[3] = { 0., 0., 0.};
 108       double interVolume = 0;
 109       for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iterTetraT = _tetra.begin(); iterTetraT != _tetra.end(); ++iterTetraT)
 110       {
 111         SplitterTetra<MyMeshType> *tmp=*iterTetraT;
 112         tmp->clearVolumesCache();
 113         double volume = tmp->intersectSourceCell(*iterCellS, baryCentre);
 114         if ( volume > 0 )
 115         {
 116           interVolume += volume;
 117           for ( int i = 0; i < SPACEDIM; ++i )
 118             total_baryCentre[i] += baryCentre[i]*volume;
 119         }
 120       }
 121       if(interVolume!=0)
 122       {
 123         for ( int i = 0; i < SPACEDIM; ++i )
 124           total_baryCentre[i] /= interVolume;
 125
 126         // coordinates of the source tetrahedron
 127         std::vector<const double*> srcCellCoords(4);
 128         for ( int n = 0; n < 4; ++n )
 129           srcCellCoords[ n ] = getCoordsOfNode( n, *iterCellS, Intersector3D<MyMeshType,MyMatrix>::_src_mesh );
 130
 131         // compute barycentric coordinates
 132         double baryCoords[4];
 133         barycentric_coords( srcCellCoords, total_baryCentre, baryCoords);
 134
 135         // store coeffs of each node of the source tetrahedron
 136         const ConnType *srcCellNodes=Intersector3D<MyMeshType,MyMatrix>::_src_mesh.getConnectivityPtr()+OTT<ConnType,numPol>::conn2C(Intersector3D<MyMeshType,MyMatrix>::_src_mesh.getConnectivityIndexPtr()[*iterCellS]);
 137         for ( int n = 0; n < 4; ++n )
 138         {
 139           double val = baryCoords[n] * interVolume;
 140           ConnType curNodeS = srcCellNodes[n];
 141           typename MyMatrix::value_type::const_iterator iterRes=resRow.find(curNodeS);
 142           if(iterRes!=resRow.end())
 143           {
 144             val += iterRes->second;
 145             resRow.erase( curNodeS );
 146           }
 147           resRow.insert(std::make_pair(curNodeS,val));
 148         }
 149       }
 150     }
 151   }
 152 }
 153
 154 #endif