1 // Copyright (C) 2007-2014 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __PlanarIntersectorP1P0Bary_TXX__
21 #define __PlanarIntersectorP1P0Bary_TXX__
23 #include "PlanarIntersectorP1P0Bary.hxx"
24 #include "InterpolationUtils.hxx"
26 #define PLAN_INTERSECTOR PlanarIntersectorP1P0Bary<MyMeshType,MyMatrix,ConcreteP1P0Intersector>
27 #define PLAN_INTER_TEMPLATE template<class MyMeshType, class MyMatrix, class ConcreteP1P0Intersector>
29 namespace INTERP_KERNEL
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)
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."
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.
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");
54 int PLAN_INTERSECTOR::getNumberOfRowsOfResMatrix() const
56 return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfElements();
60 int PLAN_INTERSECTOR::getNumberOfColsOfResMatrix() const
62 return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
66 * This method computes a value per each node of each source triangle for target.
69 void PLAN_INTERSECTOR::intersectCells(ConnType icellT,
70 const std::vector<ConnType>& icellsS,
74 std::vector<double> srcTriaCoords, tgtCellCoords, tgtCellCoordsTmp, nodeCeffs;
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();
83 typename MyMatrix::value_type& resRow=res[icellT];
85 // treat each source triangle
86 for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();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]);
93 tgtCellCoordsTmp = tgtCellCoords;
94 tgtCoords = & tgtCellCoordsTmp;
95 orientation=PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&tgtCellCoordsTmp[0], &srcTriaCoords[0],
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);
103 for(int nodeIdS=0;nodeIdS<3;nodeIdS++)
105 ConnType curNodeS=startOfCellNodeConn[nodeIdS];
106 typename MyMatrix::value_type::const_iterator iterRes=resRow.find(curNodeS);
107 if(iterRes!=resRow.end())
109 nodeCeffs[nodeIdS] += iterRes->second;
110 resRow.erase( curNodeS );
112 resRow.insert(std::make_pair(curNodeS,nodeCeffs[nodeIdS]));