From 231f7c6be4ebfd007f12bd58c066faf08ae0c057 Mon Sep 17 00:00:00 2001 From: ageay Date: Tue, 12 Nov 2013 15:11:39 +0000 Subject: [PATCH] Implementation of P0P1Bary --- src/INTERP_KERNEL/ConvexIntersector.hxx | 1 + src/INTERP_KERNEL/ConvexIntersector.txx | 1 + src/INTERP_KERNEL/Geometric2DIntersector.hxx | 1 + src/INTERP_KERNEL/Geometric2DIntersector.txx | 1 + .../PlanarIntersectorP0P1Bary.hxx | 55 +++++++++ .../PlanarIntersectorP0P1Bary.txx | 115 ++++++++++++++++++ 6 files changed, 174 insertions(+) create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.hxx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.txx diff --git a/src/INTERP_KERNEL/ConvexIntersector.hxx b/src/INTERP_KERNEL/ConvexIntersector.hxx index 3b12cdf8b..dd842b553 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.hxx +++ b/src/INTERP_KERNEL/ConvexIntersector.hxx @@ -26,6 +26,7 @@ #include "PlanarIntersectorP1P0.hxx" #include "PlanarIntersectorP1P1.hxx" #include "PlanarIntersectorP1P0Bary.hxx" +#include "PlanarIntersectorP0P1Bary.hxx" namespace INTERP_KERNEL { diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index 817677b60..f9ee8d44d 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -26,6 +26,7 @@ #include "PlanarIntersectorP1P0.txx" #include "PlanarIntersectorP1P1.txx" #include "PlanarIntersectorP1P0Bary.txx" +#include "PlanarIntersectorP0P1Bary.txx" #include "PolygonAlgorithms.txx" diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.hxx b/src/INTERP_KERNEL/Geometric2DIntersector.hxx index 5a0c9f8ff..7c077aeae 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.hxx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.hxx @@ -26,6 +26,7 @@ #include "PlanarIntersectorP1P0.hxx" #include "PlanarIntersectorP1P1.hxx" #include "PlanarIntersectorP1P0Bary.hxx" +#include "PlanarIntersectorP0P1Bary.hxx" namespace INTERP_KERNEL { diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.txx b/src/INTERP_KERNEL/Geometric2DIntersector.txx index ba4b238b1..b6bed16a9 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.txx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.txx @@ -27,6 +27,7 @@ #include "PlanarIntersectorP1P0.txx" #include "PlanarIntersectorP1P1.txx" #include "PlanarIntersectorP1P0Bary.txx" +#include "PlanarIntersectorP0P1Bary.txx" #include "CellModel.hxx" #include "InterpKernelGeo2DQuadraticPolygon.hxx" diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.hxx b/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.hxx new file mode 100644 index 000000000..31ee76537 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.hxx @@ -0,0 +1,55 @@ +// Copyright (C) 2007-2013 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 +// +// Author : Anthony Geay (CEA/DEN) + +#ifndef __PlanarIntersectorP0P1Bary_HXX__ +#define __PlanarIntersectorP0P1Bary_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP0P1Bary : public PlanarIntersector + { + public: + static const int SPACEDIM=MyMeshType::MY_SPACEDIM; + static const int MESHDIM=MyMeshType::MY_MESHDIM; + typedef typename MyMeshType::MyConnType ConnType; + static const NumberingPolicy numPol=MyMeshType::My_numPol; + protected: + PlanarIntersectorP0P1Bary(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double precision, double md3DSurf, double medianPlane, bool doRotate, int orientation, int printLevel); + public: + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + /*! + * Contrary to intersectCells method here icellS and icellT are \b not in \b C mode but in mode of MyMeshType. + */ + double intersectGeoBary(const std::vector& sourceCell, + bool sourceCellQuadratic, + const double * targetTria, + std::vector& res) + { return asLeaf().intersectGeoBary(sourceCell,sourceCellQuadratic,targetTria,res); } + protected: + ConcreteP0P1Intersector& asLeaf() { return static_cast(*this); } + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.txx b/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.txx new file mode 100644 index 000000000..b69294d2c --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1Bary.txx @@ -0,0 +1,115 @@ +// Copyright (C) 2007-2013 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 +// +// Author : Anthony Geay (CEA/DEN) + +#ifndef __PlanarIntersectorP0P1Bary_TXX__ +#define __PlanarIntersectorP0P1Bary_TXX__ + +#include "PlanarIntersectorP0P1Bary.hxx" +#include "InterpolationUtils.hxx" + +namespace INTERP_KERNEL +{ + template + PlanarIntersectorP0P1Bary::PlanarIntersectorP0P1Bary(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double precision, + double md3DSurf, double medianPlane, + bool doRotate, int orientation, int printLevel): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf, + medianPlane,doRotate,orientation,printLevel) + { + // SPEC: + // "Limitation. For the P0P1 barycentric improvement only triangle target cells in 2D and + // tetrahedrons in 3D will be supported by interpolators. If a non + // triangle/tetrahedron source cell is detected an INTERP_KERNEL::Exception should be thrown." + + // Check types of source elements here rather than in intersectCells() since a wrong type can be + // found late after a long time of calculation. + + const unsigned long numTrgElems = meshT.getNumberOfElements(); + for(unsigned long i = 0 ; i < numTrgElems ; ++i) + if ( meshT.getTypeOfElement( OTT::indFC( i )) != NORM_TRI3 ) + throw INTERP_KERNEL::Exception("P0P1 barycentric algorithm works only with triangular target meshes"); + } + + template + int PlanarIntersectorP0P1Bary::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfNodes(); + } + + template + int PlanarIntersectorP0P1Bary::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfElements(); + } + + /*! + * This method computes a value per each node of each source triangle for target. + */ + template + void PlanarIntersectorP0P1Bary::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + int orientation=1; + std::vector trgTriaCoords,trgTriaCoordsTmp; + // target cell data + PlanarIntersector::getRealTargetCoordinates(OTT::indFC(icellT),trgTriaCoords); + std::vector *tgtCoords(&trgTriaCoords); + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[icellT]); + // treat each source cells + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + std::vector srcCellCoords,srcCellCoordsTmp,nodeCeffs; + int iS=*iter; + NormalizedCellType tS=PlanarIntersector::_meshS.getTypeOfElement(OTT::indFC(iS)); + bool isSourceQuad=CellModel::GetCellModel(tS).isQuadratic(); + PlanarIntersector::getRealSourceCoordinates(OTT::indFC(iS),srcCellCoords); + std::vector *srcCoords(&srcCellCoords); + int srcNbNodes = srcCellCoords.size()/SPACEDIM; + if(SPACEDIM==3) + { + srcCellCoordsTmp=srcCellCoords; + trgTriaCoordsTmp=trgTriaCoords; + srcCoords=&srcCellCoordsTmp; + tgtCoords=&trgTriaCoordsTmp; + orientation=PlanarIntersector::projectionThis(&trgTriaCoordsTmp[0],&srcCellCoordsTmp[0], + 3,srcNbNodes); + } + //double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,targetCellCoordsTmp,isTargetQuad); + double surf=orientation*intersectGeoBary(*srcCoords,isSourceQuad,&((*tgtCoords)[0]),nodeCeffs); + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) + { + for(int nodeIdT=0;nodeIdT<3;nodeIdT++) + { + ConnType curNodeT=startOfCellNodeConn[nodeIdT]; + typename MyMatrix::value_type& resRow=res[curNodeT]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(*iter); + if(iterRes!=resRow.end()) + { + nodeCeffs[*iter] += iterRes->second; + resRow.erase(*iter); + } + resRow.insert(std::make_pair(*iter,nodeCeffs[nodeIdT])); + } + } + } + } +} +#endif -- 2.39.2