src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.txx

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