From 97863809f3466763d1269d869489bc4328f5826e Mon Sep 17 00:00:00 2001 From: ageay Date: Wed, 21 Oct 2009 06:56:27 +0000 Subject: [PATCH] Adding non conservative interpolators for 2D and 3DSurf. --- src/INTERP_KERNEL/ConvexIntersector.txx | 2 +- src/INTERP_KERNEL/InterpolationOptions.cxx | 13 +- src/INTERP_KERNEL/InterpolationOptions.hxx | 5 +- src/INTERP_KERNEL/InterpolationPlanar.txx | 41 +++- src/INTERP_KERNEL/InterpolationUtils.hxx | 7 +- src/INTERP_KERNEL/Makefile.am | 9 + src/INTERP_KERNEL/PlanarIntersectorP0P1PL.hxx | 42 ++++ src/INTERP_KERNEL/PlanarIntersectorP0P1PL.txx | 83 +++++++ src/INTERP_KERNEL/PlanarIntersectorP1P0PL.hxx | 42 ++++ src/INTERP_KERNEL/PlanarIntersectorP1P0PL.txx | 107 +++++++++ src/INTERP_KERNEL/PlanarIntersectorP1P1PL.hxx | 42 ++++ src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx | 98 ++++++++ .../PointLocator2DIntersector.hxx | 54 +++++ .../PointLocator2DIntersector.txx | 164 +++++++++++++ src/INTERP_KERNEL/PointLocatorAlgos.txx | 215 ++++++++++++++++++ src/INTERP_KERNEL/VolSurfFormulae.hxx | 14 ++ 16 files changed, 925 insertions(+), 13 deletions(-) create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP0P1PL.hxx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP0P1PL.txx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P0PL.hxx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P0PL.txx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P1PL.hxx create mode 100644 src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx create mode 100644 src/INTERP_KERNEL/PointLocator2DIntersector.hxx create mode 100644 src/INTERP_KERNEL/PointLocator2DIntersector.txx create mode 100644 src/INTERP_KERNEL/PointLocatorAlgos.txx diff --git a/src/INTERP_KERNEL/ConvexIntersector.txx b/src/INTERP_KERNEL/ConvexIntersector.txx index 317fd2393..9565e1416 100644 --- a/src/INTERP_KERNEL/ConvexIntersector.txx +++ b/src/INTERP_KERNEL/ConvexIntersector.txx @@ -177,7 +177,7 @@ namespace INTERP_KERNEL barycenter[0] = ( barycenter[0] + inter[0] + inter[SPACEDIM*(nb_inter-1)] ) / nb_inter; barycenter[1] = ( barycenter[1] + inter[1] + inter[SPACEDIM*(nb_inter-1)+1]) / nb_inter; res.resize(3); - barycentric_coords( sourceTria, &barycenter[0], &res[0]); + barycentric_coords<2>( sourceTria, &barycenter[0], &res[0]); res[0] *= area; res[1] *= area; res[2] *= area; diff --git a/src/INTERP_KERNEL/InterpolationOptions.cxx b/src/INTERP_KERNEL/InterpolationOptions.cxx index bd7cf15b0..9c5d2c31b 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.cxx +++ b/src/INTERP_KERNEL/InterpolationOptions.cxx @@ -52,6 +52,8 @@ const char INTERP_KERNEL::InterpolationOptions::CONVEX_INTERSECT2D_STR[]="Convex const char INTERP_KERNEL::InterpolationOptions::GEOMETRIC_INTERSECT2D_STR[]="Geometric2D"; +const char INTERP_KERNEL::InterpolationOptions::POINTLOCATOR_INTERSECT2D_STR[]="PointLocator2D"; + const char INTERP_KERNEL::InterpolationOptions::PLANAR_SPLIT_FACE_5_STR[]="PLANAR_FACE_5"; const char INTERP_KERNEL::InterpolationOptions::PLANAR_SPLIT_FACE_6_STR[]="PLANAR_FACE_6"; @@ -136,6 +138,11 @@ bool INTERP_KERNEL::InterpolationOptions::setOptionString(const std::string& key setIntersectionType(INTERP_KERNEL::Geometric2D); return true; } + else if(value==POINTLOCATOR_INTERSECT2D_STR) + { + setIntersectionType(INTERP_KERNEL::PointLocator2D); + return true; + } } else if(key==SPLITTING_POLICY_STR) { @@ -165,12 +172,6 @@ bool INTERP_KERNEL::InterpolationOptions::setOptionString(const std::string& key return false; } -//================================================================================ -/*! - * \brief Return "P1P0Bary" if meth=="P1P0" && this->_P1P0_bary_method=true - */ -//================================================================================ - std::string INTERP_KERNEL::InterpolationOptions::filterInterpolationMethod(const std::string& meth) const { if ( _P1P0_bary_method && meth == "P1P0" ) diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx index 2515614b1..14a5189ed 100644 --- a/src/INTERP_KERNEL/InterpolationOptions.hxx +++ b/src/INTERP_KERNEL/InterpolationOptions.hxx @@ -25,7 +25,7 @@ namespace INTERP_KERNEL { - typedef enum { Triangulation, Convex, Geometric2D } IntersectionType; + typedef enum { Triangulation, Convex, Geometric2D, PointLocator2D } IntersectionType; /// Type describing the different ways in which the hexahedron can be split into tetrahedra. /// The PLANAR_* policies persume that each face is to be considered planar, while the general /// policies make no such hypothesis. The integer at the end gives the number of tetrahedra @@ -49,7 +49,6 @@ namespace INTERP_KERNEL bool _measure_abs; SplittingPolicy _splitting_policy ; bool _P1P0_bary_method; // issue 0020440 - public: InterpolationOptions() { init(); } int getPrintLevel() const { return _print_level; } @@ -87,6 +86,7 @@ namespace INTERP_KERNEL void setP1P0BaryMethod(bool isP1P0) { _P1P0_bary_method=isP1P0; } bool getP1P0BaryMethod() const { return _P1P0_bary_method; } + std::string filterInterpolationMethod(const std::string& meth) const; void init() @@ -126,6 +126,7 @@ namespace INTERP_KERNEL static const char TRIANGULATION_INTERSECT2D_STR[]; static const char CONVEX_INTERSECT2D_STR[]; static const char GEOMETRIC_INTERSECT2D_STR[]; + static const char POINTLOCATOR_INTERSECT2D_STR[]; static const char PLANAR_SPLIT_FACE_5_STR[]; static const char PLANAR_SPLIT_FACE_6_STR[]; static const char GENERAL_SPLIT_24_STR[]; diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index 4b1729940..eca6069ae 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -30,6 +30,14 @@ #include "ConvexIntersector.txx" #include "Geometric2DIntersector.hxx" #include "Geometric2DIntersector.txx" +#include "PointLocator2DIntersector.hxx" +#include "PointLocator2DIntersector.txx" +#include "PlanarIntersectorP0P1PL.hxx" +#include "PlanarIntersectorP0P1PL.txx" +#include "PlanarIntersectorP1P0PL.hxx" +#include "PlanarIntersectorP1P0PL.txx" +#include "PlanarIntersectorP1P1PL.hxx" +#include "PlanarIntersectorP1P1PL.txx" #include "VectorUtils.hxx" #include "BBTree.txx" @@ -165,6 +173,13 @@ namespace INTERP_KERNEL InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); break; + case PointLocator2D: + intersector=new PointLocator2DIntersector(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; } } else if(meth=="P0P1") @@ -195,6 +210,13 @@ namespace INTERP_KERNEL InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); break; + case PointLocator2D: + intersector=new PlanarIntersectorP0P1PL(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; } } else if(meth=="P1P0") @@ -225,6 +247,13 @@ namespace INTERP_KERNEL InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); break; + case PointLocator2D: + intersector=new PlanarIntersectorP1P0PL(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; } } else if(meth=="P1P0Bary") @@ -255,6 +284,9 @@ namespace INTERP_KERNEL InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); break; + case PointLocator2D: + throw INTERP_KERNEL::Exception("Invalid intersector (PointLocator2D) for P1P0Bary !"); + break; } } else if(meth=="P1P1") @@ -285,10 +317,17 @@ namespace INTERP_KERNEL InterpolationOptions::getPrecision(), InterpolationOptions::getOrientation()); break; + case PointLocator2D: + intersector=new PlanarIntersectorP1P1PL(myMeshT, myMeshS, _dim_caracteristic, + InterpolationOptions::getMaxDistance3DSurfIntersect(), + InterpolationOptions::getMedianPlane(), + InterpolationOptions::getPrecision(), + InterpolationOptions::getOrientation()); + break; } } else - throw INTERP_KERNEL::Exception("Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\""); + throw INTERP_KERNEL::Exception("Invalid method specified or intersection type ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\""); /****************************************************************/ /* Create a search tree based on the bounding boxes */ /* Instanciate the intersector and initialise the result vector */ diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index 3198b5cc2..ed49f0e4b 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -213,12 +213,13 @@ namespace INTERP_KERNEL /* triaCoords are in full interlace */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ + template inline void barycentric_coords(const double* triaCoords, const double* p, double* bc) { // matrix 2x2 double - T11 = triaCoords[0]-triaCoords[4], T12 = triaCoords[2]-triaCoords[4], - T21 = triaCoords[1]-triaCoords[5], T22 = triaCoords[3]-triaCoords[5]; + T11 = triaCoords[0]-triaCoords[2*SPACEDIM], T12 = triaCoords[SPACEDIM]-triaCoords[2*SPACEDIM], + T21 = triaCoords[1]-triaCoords[2*SPACEDIM+1], T22 = triaCoords[SPACEDIM+1]-triaCoords[2*SPACEDIM+1]; // matrix determinant double Tdet = T11*T22 - T12*T21; if ( fabs( Tdet ) < std::numeric_limits::min() ) { @@ -228,7 +229,7 @@ namespace INTERP_KERNEL // matrix inverse double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11; // vector - double r11 = p[0]-triaCoords[4], r12 = p[1]-triaCoords[5]; + double r11 = p[0]-triaCoords[2*SPACEDIM], r12 = p[1]-triaCoords[2*SPACEDIM+1]; // barycentric coordinates: mutiply matrix by vector bc[0] = (t11 * r11 + t12 * r12)/Tdet; bc[1] = (t21 * r11 + t22 * r12)/Tdet; diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 4cba364f1..ca8391823 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -38,6 +38,8 @@ ConvexIntersector.hxx \ ConvexIntersector.txx \ Geometric2DIntersector.hxx \ Geometric2DIntersector.txx \ +PointLocator2DIntersector.hxx \ +PointLocator2DIntersector.txx \ INTERPKERNELDefines.hxx \ InterpKernelMatrix.hxx \ Interpolation.hxx \ @@ -72,6 +74,7 @@ MeshElement.txx \ MeshRegion.hxx \ MeshRegion.txx \ MeshUtils.hxx \ +PointLocatorAlgos.txx \ PlanarIntersector.hxx \ PlanarIntersector.txx \ PlanarIntersectorP0P0.hxx \ @@ -82,6 +85,12 @@ PlanarIntersectorP1P0.hxx \ PlanarIntersectorP1P0.txx \ PlanarIntersectorP1P0Bary.hxx \ PlanarIntersectorP1P0Bary.txx \ +PlanarIntersectorP0P1PL.hxx \ +PlanarIntersectorP0P1PL.txx \ +PlanarIntersectorP1P0PL.hxx \ +PlanarIntersectorP1P0PL.txx \ +PlanarIntersectorP1P1PL.hxx \ +PlanarIntersectorP1P1PL.txx \ PlanarIntersectorP1P1.hxx \ PlanarIntersectorP1P1.txx \ PolygonAlgorithms.hxx \ diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.hxx b/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.hxx new file mode 100644 index 000000000..ce39816ed --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.hxx @@ -0,0 +1,42 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP0P1PL_HXX__ +#define __PLANARINTERSECTORP0P1PL_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP0P1PL : 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; + public: + PlanarIntersectorP0P1PL(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation); + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.txx b/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.txx new file mode 100644 index 000000000..2f71ccc85 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP0P1PL.txx @@ -0,0 +1,83 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP0P1PL_TXX__ +#define __PLANARINTERSECTORP0P1PL_TXX__ + +#include "PlanarIntersectorP0P1PL.hxx" +#include "PlanarIntersector.txx" +#include "CellModel.hxx" + +#include "PointLocatorAlgos.txx" + +namespace INTERP_KERNEL +{ + template + PlanarIntersectorP0P1PL::PlanarIntersectorP0P1PL(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, + double medianPlane, double precision, int orientation): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,true,orientation,0) + { + } + + template + void PlanarIntersectorP0P1PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector< std::vector > coordsOfSources(icellsS.size()); + int ii=0; + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++,ii++) + PlanarIntersector::getRealSourceCoordinates(*iter,coordsOfSources[ii]); + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[icellT]); + std::vector coordsTarget; + PlanarIntersector::getRealTargetCoordinates(icellT,coordsTarget); + int nbNodesT=coordsTarget.size()/SPACEDIM; + ii=0; + for(typename std::vector::const_iterator iter2=icellsS.begin();iter2!=icellsS.end();iter2++,ii++) + { + std::vector tmpSource(coordsOfSources[ii]); + std::vector tmpTarget(coordsTarget); + if(SPACEDIM==3) + PlanarIntersector::projectionThis(&tmpSource[0],&tmpTarget[0],tmpSource.size()/SPACEDIM,nbNodesT); + for(int nodeIdT=0;nodeIdT::isElementContainsPointAlg2D(&tmpTarget[0]+nodeIdT*SPACEDIM,&tmpSource[0],tmpSource.size()/SPACEDIM) ) + { + ConnType curNodeTInCmode=OTT::coo2C(startOfCellNodeConn[nodeIdT]); + typename MyMatrix::value_type& resRow=res[curNodeTInCmode]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(*iter2)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(*iter2),1.)); + } + } + } + } + + template + int PlanarIntersectorP0P1PL::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfNodes(); + } + + template + int PlanarIntersectorP0P1PL::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfElements(); + } +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.hxx new file mode 100644 index 000000000..dd78a7803 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.hxx @@ -0,0 +1,42 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP1P0PL_HXX__ +#define __PLANARINTERSECTORP1P0PL_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP1P0PL : 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; + public: + PlanarIntersectorP1P0PL(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation); + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.txx new file mode 100644 index 000000000..fbc655c9c --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0PL.txx @@ -0,0 +1,107 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP1P0PL_TXX__ +#define __PLANARINTERSECTORP1P0PL_TXX__ + +#include "PlanarIntersectorP1P0PL.hxx" +#include "PlanarIntersector.txx" +#include "CellModel.hxx" + +#include "PointLocatorAlgos.txx" +#include "QuadraticPolygon.hxx" +#include "MeshUtils.hxx" + +namespace INTERP_KERNEL +{ + template + PlanarIntersectorP1P0PL::PlanarIntersectorP1P0PL(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, + double medianPlane, double precision, int orientation): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,true,orientation,0) + { + } + + template + void PlanarIntersectorP1P0PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector CoordsT; + typename MyMatrix::value_type& resRow=res[icellT]; + PlanarIntersector::getRealTargetCoordinates(icellT,CoordsT); + double baryT[SPACEDIM]; + double baryTTmp[SPACEDIM]; + calculateBarycenterDyn2(&CoordsT[0],CoordsT.size()/SPACEDIM,baryT); + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + NormalizedCellType tS=PlanarIntersector::_meshS.getTypeOfElement(*iter); + if(tS!=NORM_TRI3) + throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==2. Only TRI3 supported !"); + std::vector CoordsS; + PlanarIntersector::getRealSourceCoordinates(*iter,CoordsS); + if(SPACEDIM==2) + { + std::copy(baryT,baryT+SPACEDIM,baryTTmp); + } + else + { + double littleTargetCell[9]; + std::copy(baryT,baryT+3,littleTargetCell); + std::copy(CoordsT.begin(),CoordsT.begin()+3,littleTargetCell+3); + std::copy(CoordsT.begin()+3,CoordsT.begin()+6,littleTargetCell+6); + PlanarIntersector::projectionThis(&CoordsS[0],littleTargetCell,3,3); + std::copy(littleTargetCell,littleTargetCell+3,baryTTmp); + } + if( PointLocatorAlgos::isElementContainsPointAlg2D(baryTTmp,&CoordsS[0],3) ) + { + double resLoc[3]; + barycentric_coords(&CoordsS[0],baryTTmp,resLoc); + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[*iter]); + for(int nodeIdS=0;nodeIdS<3;nodeIdS++) + { + if(fabs(resLoc[nodeIdS])>PlanarIntersector::_precision) + { + ConnType curNodeSInCmode=OTT::coo2C(startOfCellNodeConn[nodeIdS]); + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),resLoc[nodeIdS])); + else + { + double val=(*iterRes).second+resLoc[nodeIdS]; + resRow.erase(OTT::indFC(curNodeSInCmode)); + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),val)); + } + } + } + } + } + } + + template + int PlanarIntersectorP1P0PL::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfElements(); + } + + template + int PlanarIntersectorP1P0PL::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfNodes(); + } +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.hxx new file mode 100644 index 000000000..8f82b5495 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.hxx @@ -0,0 +1,42 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP1P1PL_HXX__ +#define __PLANARINTERSECTORP1P1PL_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP1P1PL : 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; + public: + PlanarIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation); + void intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx new file mode 100644 index 000000000..294091f91 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P1PL.txx @@ -0,0 +1,98 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __PLANARINTERSECTORP1P1PL_TXX__ +#define __PLANARINTERSECTORP1P1PL_TXX__ + +#include "PlanarIntersectorP1P1PL.hxx" +#include "PlanarIntersector.txx" +#include "CellModel.hxx" + +#include "PointLocatorAlgos.txx" +#include "MeshUtils.hxx" + +namespace INTERP_KERNEL +{ + template + PlanarIntersectorP1P1PL::PlanarIntersectorP1P1PL(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, + double medianPlane, double precision, int orientation): + PlanarIntersector(meshT,meshS,dimCaracteristic,precision,md3DSurf,medianPlane,true,orientation,0) + { + } + + template + void PlanarIntersectorP1P1PL::intersectCells(ConnType icellT, const std::vector& icellsS, MyMatrix& res) + { + std::vector CoordsT; + PlanarIntersector::getRealTargetCoordinates(icellT,CoordsT); + int nbOfNodesT=CoordsT.size()/SPACEDIM; + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + NormalizedCellType tS=PlanarIntersector::_meshS.getTypeOfElement(*iter); + if(tS!=NORM_TRI3) + throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==2. Only TRI3 supported !"); + std::vector CoordsS; + PlanarIntersector::getRealSourceCoordinates(*iter,CoordsS); + std::vector CoordsTTmp(CoordsT); + if(SPACEDIM==3) + PlanarIntersector::projectionThis(&CoordsS[0],&CoordsTTmp[0],CoordsS.size()/SPACEDIM,nbOfNodesT); + const ConnType *startOfCellNodeConnT=PlanarIntersector::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[icellT]); + for(int nodeIdT=0;nodeIdT::ind2C(startOfCellNodeConnT[nodeIdT])]; + if( PointLocatorAlgos::isElementContainsPointAlg2D(&CoordsTTmp[nodeIdT*SPACEDIM],&CoordsS[0],3) ) + { + double resLoc[3]; + barycentric_coords(&CoordsS[0],&CoordsTTmp[nodeIdT*SPACEDIM],resLoc); + const ConnType *startOfCellNodeConnS=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[*iter]); + for(int nodeIdS=0;nodeIdS<3;nodeIdS++) + { + if(fabs(resLoc[nodeIdS])>PlanarIntersector::_precision) + { + ConnType curNodeSInCmode=OTT::coo2C(startOfCellNodeConnS[nodeIdS]); + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT::indFC(curNodeSInCmode)); + if(iterRes==resRow.end()) + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),resLoc[nodeIdS])); + else + { + double val=(*iterRes).second+resLoc[nodeIdS]; + resRow.erase(OTT::indFC(curNodeSInCmode)); + resRow.insert(std::make_pair(OTT::indFC(curNodeSInCmode),val)); + } + } + } + } + } + } + } + + template + int PlanarIntersectorP1P1PL::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfNodes(); + } + + template + int PlanarIntersectorP1P1PL::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfNodes(); + } +} + +#endif diff --git a/src/INTERP_KERNEL/PointLocator2DIntersector.hxx b/src/INTERP_KERNEL/PointLocator2DIntersector.hxx new file mode 100644 index 000000000..e41841ba2 --- /dev/null +++ b/src/INTERP_KERNEL/PointLocator2DIntersector.hxx @@ -0,0 +1,54 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __POINTLOCATORINTERSECTOR_HXX__ +#define __POINTLOCATORINTERSECTOR_HXX__ + +#include "PlanarIntersectorP0P0.hxx" +#include "PlanarIntersectorP0P1.hxx" +#include "PlanarIntersectorP1P0.hxx" +#include "PlanarIntersectorP1P1.hxx" +#include "PlanarIntersectorP1P0Bary.hxx" + +namespace INTERP_KERNEL +{ + class QuadraticPolygon; + + template class InterpType> + class PointLocator2DIntersector : public InterpType > + { + 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; + public: + PointLocator2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, double medianPlane, double precision, int orientation); + double intersectGeometry(ConnType icellT, ConnType icellS, ConnType nbNodesT, ConnType nbNodesS); + double intersectGeometryWithQuadrangle(const double *quadrangle, const std::vector& sourceCoords, bool isSourceQuad); + double intersectGeometryGeneral(const std::vector& targetCoords, const std::vector& sourceCoords); + double intersectGeoBary(const std::vector& targetCell, bool targetCellQuadratic, const double *sourceCell, std::vector& res); + private: + QuadraticPolygon *buildPolygonFrom(const std::vector& coords, NormalizedCellType type); + QuadraticPolygon *buildPolygonAFrom(ConnType cell, int nbOfPoints, NormalizedCellType type); + QuadraticPolygon *buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type); + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PointLocator2DIntersector.txx b/src/INTERP_KERNEL/PointLocator2DIntersector.txx new file mode 100644 index 000000000..c88e9be52 --- /dev/null +++ b/src/INTERP_KERNEL/PointLocator2DIntersector.txx @@ -0,0 +1,164 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __POINTLOCATORINTERSECTOR_TXX__ +#define __POINTLOCATORINTERSECTOR_TXX__ + +#include "PointLocator2DIntersector.hxx" +#include "PlanarIntersectorP0P0.txx" +#include "PlanarIntersectorP0P1.txx" +#include "PlanarIntersectorP1P0.txx" +#include "PlanarIntersectorP1P1.txx" +#include "PlanarIntersectorP1P0Bary.txx" +#include "CellModel.hxx" + +#include "QuadraticPolygon.hxx" +#include "PointLocatorAlgos.txx" + +#define PTLOC2D_INTERSECTOR PointLocator2DIntersector +#define INTERSECTOR_TEMPLATE template class InterpType> + +namespace INTERP_KERNEL +{ + INTERSECTOR_TEMPLATE + PTLOC2D_INTERSECTOR::PointLocator2DIntersector(const MyMeshType& meshT, const MyMeshType& meshS, + double dimCaracteristic, double md3DSurf, double medianPlane, + double precision, int orientation): + InterpType(meshT,meshS,dimCaracteristic, precision, md3DSurf, medianPlane, true, orientation, 0) + { + } + + INTERSECTOR_TEMPLATE + double PTLOC2D_INTERSECTOR::intersectGeometry(ConnType icellT, ConnType icellS, + ConnType nbNodesT, ConnType nbNodesS) + { + int orientation = 1; + std::vector CoordsT; + std::vector CoordsS; + PlanarIntersector::getRealCoordinates(icellT,icellS,nbNodesT,nbNodesS,CoordsT,CoordsS,orientation); + NormalizedCellType tT=PlanarIntersector::_meshT.getTypeOfElement(icellT); + QuadraticPolygon *pT=buildPolygonFrom(CoordsT,tT); + double baryT[SPACEDIM]; + pT->getBarycenter(baryT); + delete pT; + if( PointLocatorAlgos::isElementContainsPointAlg2D(baryT,&CoordsS[0],nbNodesS) ) + return 1.; + return 0.; + } + + INTERSECTOR_TEMPLATE + double PTLOC2D_INTERSECTOR::intersectGeometryWithQuadrangle(const double * quadrangle, + const std::vector& sourceCoords, + bool isSourceQuad) + { + int nbOfSourceNodes=sourceCoords.size()/SPACEDIM; + std::vector nodes2(nbOfSourceNodes); + for(int i=0;igetBarycenter(bary); + delete p2; + if( PointLocatorAlgos::isElementContainsPointAlg2D(bary,quadrangle,4) ) + return 1.; + return 0.; + } + + INTERSECTOR_TEMPLATE + double PTLOC2D_INTERSECTOR::intersectGeometryGeneral(const std::vector& targetCoords, + const std::vector& sourceCoords) + { + int nbOfTargetNodes=targetCoords.size()/SPACEDIM; + int nbOfSourceNodes=sourceCoords.size()/SPACEDIM; + std::vector nodes2(nbOfSourceNodes); + for(int i=0;igetBarycenter(bary); + delete p; + if( PointLocatorAlgos::isElementContainsPointAlg2D(bary,&targetCoords[0],nbOfTargetNodes) ) + return 1.; + return 0.; + } + + //================================================================================ + /*! + * \brief Intersect a triangle and a polygon for P1P0 barycentric algorithm + * \param targetCell - list of coordinates of target polygon in full interlace + * \param targetCellQuadratic - specifies if target polygon is quadratic or not + * \param sourceTria - list of coordinates of source triangle + * \param res - coefficients a,b and c associated to nodes of sourceTria + */ + //================================================================================ + + INTERSECTOR_TEMPLATE + double PTLOC2D_INTERSECTOR::intersectGeoBary(const std::vector& targetCell, + bool targetCellQuadratic, + const double * sourceTria, + std::vector& res) + { + throw INTERP_KERNEL::Exception("intersectGeoBary incompatible with PointLocator. Desactivate P1P0Bary to avoid the problem"); + return 0.; + } + + INTERSECTOR_TEMPLATE + QuadraticPolygon *PTLOC2D_INTERSECTOR::buildPolygonFrom(const std::vector& coords, NormalizedCellType type) + { + int nbNodes=coords.size()/SPACEDIM; + std::vector nodes(nbNodes); + for(int i=0;i::_connectT+OTT::conn2C(PlanarIntersector::_connIndexT[OTT::ind2C(cell)]); + std::vector nodes(nbOfPoints); + for(int i=0;i::_coordsT+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); + if(CellModel::getCellModel(type).isQuadratic()) + return QuadraticPolygon::buildLinearPolygon(nodes); + else + return QuadraticPolygon::buildArcCirclePolygon(nodes); + } + + INTERSECTOR_TEMPLATE + QuadraticPolygon *PTLOC2D_INTERSECTOR::buildPolygonBFrom(ConnType cell, int nbOfPoints, NormalizedCellType type) + { + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[OTT::ind2C(cell)]); + std::vector nodes(nbOfPoints); + for(int i=0;i::_coordsS+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); + if(type!=NORM_TRI6 && type!=NORM_QUAD8) + return QuadraticPolygon::buildLinearPolygon(nodes); + else + return QuadraticPolygon::buildArcCirclePolygon(nodes); + } +} + +#endif diff --git a/src/INTERP_KERNEL/PointLocatorAlgos.txx b/src/INTERP_KERNEL/PointLocatorAlgos.txx new file mode 100644 index 000000000..c4b4cdf35 --- /dev/null +++ b/src/INTERP_KERNEL/PointLocatorAlgos.txx @@ -0,0 +1,215 @@ +// Copyright (C) 2007-2008 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 +// +#ifndef __POINTLOCATORALGOS_TXX__ +#define __POINTLOCATORALGOS_TXX__ + +#include "InterpolationUtils.hxx" +#include "CellModel.hxx" +#include "BBTree.txx" + +#include +#include + +namespace INTERP_KERNEL +{ + class GenericPointLocatorAlgos + { + public: + virtual ~GenericPointLocatorAlgos() { } + virtual std::list locates(const double* x) = 0; + }; + + template + class PointLocatorAlgos: public GenericPointLocatorAlgos + { + private : + double* _bb; + BBTree* _tree; + const MyMeshType& _mesh; + public: + PointLocatorAlgos(const MyMeshType& mesh):_mesh(mesh) + { + typedef typename MyMeshType::MyConnType ConnType; + const int SPACEDIM=MyMeshType::MY_SPACEDIM; + const NumberingPolicy numPol=MyMeshType::My_numPol; + int nelem = _mesh.getNumberOfElements(); + _bb = new double[SPACEDIM*2*nelem]; + const ConnType* conn = _mesh.getConnectivityPtr(); + const ConnType* conn_index = _mesh.getConnectivityIndexPtr(); + const double* coords=_mesh.getCoordinatesPtr(); + for (int i=0; i::max(); + _bb[2*(i*SPACEDIM+idim)+1]=-std::numeric_limits::max(); + } + for (int index= conn_index[i]; index < conn_index[i+1];index++) + { + //coordelem points to the coordinates of the current node of the i-th element + const double* coordelem = coords+OTT::ind2C(conn[OTT::ind2C(index)])*SPACEDIM; + + //the bounding box is updated by checking wheher the node is at the min/max in exach dimension + for (int idim=0; idim_bb[2*(i*SPACEDIM+idim)+1])?coordelem[idim]:_bb[2*(i*SPACEDIM+idim)+1]; + } + } + } + _tree=new BBTree(_bb,0,0,nelem); + } + + ~PointLocatorAlgos() + { + delete[] _bb; + delete _tree; + } + + //returns the list of elements that contains + //the point pointed to by x + std::list locates(const double* x) + { + typedef typename MyMeshType::MyConnType ConnType; + const NumberingPolicy numPol=MyMeshType::My_numPol; + std::vector candidates; + _tree->getElementsAroundPoint(x,candidates); + std::list retlist; + for(int i=0; i< candidates.size(); i++) + { + int ielem=candidates[i]; + if (elementContainsPoint(ielem,x)) + retlist.push_back(OTT::indFC(ielem)); + } + return retlist; + } + + static bool isElementContainsPointAlg2D(const double *ptToTest, const double *cellPts, int nbEdges) + { + // with dimension 2, it suffices to check all the edges + // and see if the sign of double products from the point + //is always the same. + // C + // / \ + // / \ + // Xo / \ + // A-------B + // + //here XA^XC and XC^XB have different signs + const int SPACEDIM=MyMeshType::MY_SPACEDIM; + int* sign = new int[nbEdges]; + for (int iedge=0; iedge1e-12) + sign[iedge]=1; + else + sign[iedge]=0; + } + bool ret=decideFromSign(sign, nbEdges); + delete [] sign; + return ret; + } + + static bool isElementContainsPointAlg3D(const double *ptToTest, const double *cellPts, int nbEdges) + { + return true; + } + + static bool isElementContainsPoint(const double *ptToTest, NormalizedCellType type, const double *coords, const int *conn_elem) + { + const int SPACEDIM=MyMeshType::MY_SPACEDIM; + typedef typename MyMeshType::MyConnType ConnType; + const NumberingPolicy numPol=MyMeshType::My_numPol; + + const CellModel& cmType=CellModel::getCellModel(type); + // + if (SPACEDIM==2) + { + int nbEdges=cmType.getNumberOfSons(); + double *pts = new double[nbEdges*SPACEDIM]; + for (int iedge=0; iedge::ind2C(conn_elem[iedge])); + std::copy(a,a+SPACEDIM,pts+iedge*SPACEDIM); + } + bool ret=isElementContainsPointAlg2D(ptToTest,pts,nbEdges); + delete [] pts; + return ret; + } + + if (SPACEDIM==3) + { + int nbfaces = cmType.getNumberOfSons(); + int* sign = new int[nbfaces]; + for (int iface=0; iface::ind2C(conn_elem[connface[0]])); + const double* BB=coords+SPACEDIM*(OTT::ind2C(conn_elem[connface[1]])); + const double* CC=coords+SPACEDIM*(OTT::ind2C(conn_elem[connface[2]])); + + double Vol=triple_product(AA,BB,CC,ptToTest); + if (Vol<-1e-12) + sign[iface]=-1; + else if (Vol>1e-12) + sign[iface]=1; + else + sign[iface]=0; + } + bool ret=decideFromSign(sign, nbfaces); + delete [] sign; + return ret; + } + } + + bool elementContainsPoint(typename MyMeshType::MyConnType i, const double* x) + { + //as i is extracted from the BBTRee, it is already in C numbering + //it is not necessary to convert it from F to C + typedef typename MyMeshType::MyConnType ConnType; + const NumberingPolicy numPol=MyMeshType::My_numPol; + + const double* coords= _mesh.getCoordinatesPtr(); + const ConnType* conn=_mesh.getConnectivityPtr(); + const ConnType* conn_index= _mesh.getConnectivityIndexPtr(); + const ConnType* conn_elem=conn+OTT::ind2C(conn_index[i]); + NormalizedCellType type=_mesh.getTypeOfElement(OTT::indFC(i)); + return isElementContainsPoint(x,type,coords,conn_elem); + } + + static bool decideFromSign(const int* sign, int nbelem) + { + int min_sign = 1; + int max_sign = -1; + for (int i=0; imax_sign)?sign[i]:max_sign; + } + return (min_sign!=-1 || max_sign!=1); + } + }; +} + +#endif diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx index 6e4a3d18a..74b4cf32d 100644 --- a/src/INTERP_KERNEL/VolSurfFormulae.hxx +++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx @@ -488,6 +488,20 @@ namespace INTERP_KERNEL bary[i]=temp/nbPts; } } + + template + inline void calculateBarycenterDyn2(const double *pts, int nbPts, double *bary) + { + for(int i=0;i