src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.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 __PlanarIntersectorP1P0Bary_TXX__
  21 #define __PlanarIntersectorP1P0Bary_TXX__
  22
  23 #include "PlanarIntersectorP1P0Bary.hxx"
  24 #include "InterpolationUtils.hxx"
  25
  26 #define PLAN_INTERSECTOR PlanarIntersectorP1P0Bary<MyMeshType,MyMatrix,ConcreteP1P0Intersector>
  27 #define PLAN_INTER_TEMPLATE template<class MyMeshType, class MyMatrix, class ConcreteP1P0Intersector>
  28
  29 namespace INTERP_KERNEL
  30 {
  31   PLAN_INTER_TEMPLATE
  32   PLAN_INTERSECTOR::PlanarIntersectorP1P0Bary(const MyMeshType& meshT, const MyMeshType& meshS,
  33                                               double dimCaracteristic, double precision,
  34                                               double md3DSurf, double minDot3DSurf, double medianPlane,
  35                                               bool doRotate, int orientation, int printLevel):
  36     PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,minDot3DSurf,
  37                                            medianPlane,doRotate,orientation,printLevel)
  38   {
  39     // SPEC:
  40     // "Limitation. For the P1P0 barycentric improvement only triangle source cells in 2D and
  41     // tetrahedrons in 3D will be supported by interpolators. If a non
  42     // triangle/tetrahedron source cell is detected an INTERP_KERNEL::Exception should be thrown."
  43
  44     // Check types of source elements here rather than in intersectCells() since a wrong type can be
  45     // found late after a long time of calculation.
  46
  47     const unsigned long numSrcElems = meshS.getNumberOfElements();
  48     for(unsigned long i = 0 ; i < numSrcElems ; ++i)
  49       if ( meshS.getTypeOfElement( OTT<ConnType,numPol>::indFC( i )) != NORM_TRI3 )
  50         throw INTERP_KERNEL::Exception("P1P0 barycentric algorithm works only with triangular source meshes");
  51   }
  52
  53   PLAN_INTER_TEMPLATE
  54   int PLAN_INTERSECTOR::getNumberOfRowsOfResMatrix() const
  55   {
  56     return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfElements();
  57   }
  58
  59   PLAN_INTER_TEMPLATE
  60   int PLAN_INTERSECTOR::getNumberOfColsOfResMatrix() const
  61   {
  62     return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
  63   }
  64
  65   /*!
  66    * This method computes a value per each node of each source triangle for target.
  67    */
  68   PLAN_INTER_TEMPLATE
  69   void PLAN_INTERSECTOR::intersectCells(ConnType                     icellT,
  70                                         const std::vector<ConnType>& icellsS,
  71                                         MyMatrix&                    res)
  72   {
  73     int orientation=1;
  74     std::vector<double> srcTriaCoords, tgtCellCoords, tgtCellCoordsTmp, nodeCeffs;
  75
  76     // target cell data
  77     PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),tgtCellCoords);
  78     std::vector<double> * tgtCoords = & tgtCellCoords;
  79     int tgtNbNodes = tgtCellCoords.size()/SPACEDIM;
  80     NormalizedCellType tT=PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getTypeOfElement(OTT<ConnType,numPol>::indFC(icellT));
  81     bool isTargetQuad=CellModel::GetCellModel(tT).isQuadratic();
  82
  83     typename MyMatrix::value_type& resRow=res[icellT];
  84
  85     // treat each source triangle
  86     for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
  87     {
  88       int iS=*iter;
  89       PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(OTT<ConnType,numPol>::indFC(iS),srcTriaCoords);
  90       const ConnType *startOfCellNodeConn=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[iS]);
  91       if(SPACEDIM==3)
  92       {
  93         tgtCellCoordsTmp = tgtCellCoords;
  94         tgtCoords = & tgtCellCoordsTmp;
  95         orientation=PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&tgtCellCoordsTmp[0], &srcTriaCoords[0],
  96                                                                            tgtNbNodes, 3);
  97       }
  98       //double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,targetCellCoordsTmp,isTargetQuad);
  99       double surf=orientation*intersectGeoBary( *tgtCoords, isTargetQuad, &srcTriaCoords[0], nodeCeffs );
 100       surf=PlanarIntersector<MyMeshType,MyMatrix>::getValueRegardingOption(surf);
 101       if(surf!=0.)
 102       {
 103         for(int nodeIdS=0;nodeIdS<3;nodeIdS++)
 104         {
 105           ConnType curNodeS=startOfCellNodeConn[nodeIdS];
 106           typename MyMatrix::value_type::const_iterator iterRes=resRow.find(curNodeS);
 107           if(iterRes!=resRow.end())
 108           {
 109             nodeCeffs[nodeIdS] += iterRes->second;
 110             resRow.erase( curNodeS );
 111           }
 112           resRow.insert(std::make_pair(curNodeS,nodeCeffs[nodeIdS]));
 113         }
 114       }
 115     }
 116   }
 117 }
 118 #endif