From: eap Date: Tue, 29 Sep 2009 06:48:54 +0000 (+0000) Subject: 0020440: [CEA 349] P1P0 barycentric interpolators X-Git-Tag: V5_1_main_FINAL~343 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=06271937a8e2beeedbbbcf3a959558f49482f577;p=tools%2Fmedcoupling.git 0020440: [CEA 349] P1P0 barycentric interpolators +Intersector3DP1P0Bary.hxx \ +Intersector3DP1P0Bary.txx \ +PlanarIntersectorP1P0Bary.hxx \ +PlanarIntersectorP1P0Bary.txx \ +PolyhedronIntersectorP1P0Bary.hxx \ +PolyhedronIntersectorP1P0Bary.txx \ --- diff --git a/src/INTERP_KERNEL/Intersector3DP1P0Bary.hxx b/src/INTERP_KERNEL/Intersector3DP1P0Bary.hxx new file mode 100644 index 000000000..c2eb30314 --- /dev/null +++ b/src/INTERP_KERNEL/Intersector3DP1P0Bary.hxx @@ -0,0 +1,36 @@ +// 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 __Intersector3DP1P0Bary_HXX__ +#define __Intersector3DP1P0Bary_HXX__ + +#include "Intersector3D.hxx" + +namespace INTERP_KERNEL +{ + template + class Intersector3DP1P0Bary : public Intersector3D + { + public: + Intersector3DP1P0Bary(const MyMeshType& targetMesh, const MyMeshType& srcMesh); + int getNumberOfRowsOfResMatrix() const; + int getNumberOfColsOfResMatrix() const; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Intersector3DP1P0Bary.txx b/src/INTERP_KERNEL/Intersector3DP1P0Bary.txx new file mode 100644 index 000000000..9e94aca9c --- /dev/null +++ b/src/INTERP_KERNEL/Intersector3DP1P0Bary.txx @@ -0,0 +1,45 @@ +// 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 __Intersector3DP1P0Bary_TXX__ +#define __Intersector3DP1P0Bary_TXX__ + +#include "Intersector3DP1P0Bary.hxx" +#include "Intersector3D.txx" + +namespace INTERP_KERNEL +{ + template + Intersector3DP1P0Bary::Intersector3DP1P0Bary(const MyMeshType& targetMesh, const MyMeshType& srcMesh):Intersector3D(targetMesh,srcMesh) + { + } + + template + int Intersector3DP1P0Bary::getNumberOfRowsOfResMatrix() const + { + return Intersector3D::_target_mesh.getNumberOfElements(); + } + + template + int Intersector3DP1P0Bary::getNumberOfColsOfResMatrix() const + { + return Intersector3D::_src_mesh.getNumberOfNodes(); + } +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 3ca8017a9..1492b0730 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -60,6 +60,8 @@ Intersector3DP0P1.hxx \ Intersector3DP0P1.txx \ Intersector3DP1P0.hxx \ Intersector3DP1P0.txx \ +Intersector3DP1P0Bary.hxx \ +Intersector3DP1P0Bary.txx \ Log.hxx \ MeshElement.hxx \ MeshElement.txx \ @@ -74,6 +76,8 @@ PlanarIntersectorP0P1.hxx \ PlanarIntersectorP0P1.txx \ PlanarIntersectorP1P0.hxx \ PlanarIntersectorP1P0.txx \ +PlanarIntersectorP1P0Bary.hxx \ +PlanarIntersectorP1P0Bary.txx \ PlanarIntersectorP1P1.hxx \ PlanarIntersectorP1P1.txx \ PolygonAlgorithms.hxx \ @@ -84,6 +88,8 @@ PolyhedronIntersectorP0P1.hxx \ PolyhedronIntersectorP0P1.txx \ PolyhedronIntersectorP1P0.hxx \ PolyhedronIntersectorP1P0.txx \ +PolyhedronIntersectorP1P0Bary.hxx \ +PolyhedronIntersectorP1P0Bary.txx \ RegionNode.hxx \ SplitterTetra.hxx \ SplitterTetra.txx \ @@ -101,6 +107,11 @@ VTKNormalizedUnstructuredMesh.hxx \ VTKNormalizedUnstructuredMesh.txx \ VectorUtils.hxx +#PolyhedronIntersectorP1P1.hxx +#PolyhedronIntersectorP1P1.txx +#Intersector3DP1P1.hxx +#Intersector3DP1P1.txx + EXTRA_DIST += \ InterpKernelUtilities.hxx \ Intersector3DP0P0.hxx \ diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.hxx b/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.hxx new file mode 100644 index 000000000..8a8a95702 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.hxx @@ -0,0 +1,53 @@ +// 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 __PlanarIntersectorP1P0Bary_HXX__ +#define __PlanarIntersectorP1P0Bary_HXX__ + +#include "PlanarIntersector.hxx" + +namespace INTERP_KERNEL +{ + template + class PlanarIntersectorP1P0Bary : 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: + PlanarIntersectorP1P0Bary(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& targetCell, + bool targetCellQuadratic, + const double * sourceTria, + std::vector& res) + { return asLeaf().intersectGeoBary(targetCell,targetCellQuadratic,sourceTria,res); } + protected: + ConcreteP1P0Intersector& asLeaf() { return static_cast(*this); } + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.txx b/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.txx new file mode 100644 index 000000000..6e6d666f0 --- /dev/null +++ b/src/INTERP_KERNEL/PlanarIntersectorP1P0Bary.txx @@ -0,0 +1,104 @@ +// 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 __PlanarIntersectorP1P0Bary_TXX__ +#define __PlanarIntersectorP1P0Bary_TXX__ + +#include "PlanarIntersectorP1P0Bary.hxx" +#include "InterpolationUtils.hxx" + +#define PLAN_INTERSECTOR PlanarIntersectorP1P0Bary +#define PLAN_INTER_TEMPLATE template + +namespace INTERP_KERNEL +{ + PLAN_INTER_TEMPLATE + PLAN_INTERSECTOR::PlanarIntersectorP1P0Bary(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) + { + } + + PLAN_INTER_TEMPLATE + int PLAN_INTERSECTOR::getNumberOfRowsOfResMatrix() const + { + return PlanarIntersector::_meshT.getNumberOfElements(); + } + + PLAN_INTER_TEMPLATE + int PLAN_INTERSECTOR::getNumberOfColsOfResMatrix() const + { + return PlanarIntersector::_meshS.getNumberOfNodes(); + } + + /*! + * This method computes a value per each node of each source triangle for target. + */ + PLAN_INTER_TEMPLATE + void PLAN_INTERSECTOR::intersectCells(ConnType icellT, + const std::vector& icellsS, + MyMatrix& res) + { + int orientation=1; + std::vector srcTriaCoords, tgtCellCoords, tgtCellCoordsTmp, nodeCeffs; + + // target cell data + PlanarIntersector::getRealTargetCoordinates(OTT::indFC(icellT),tgtCellCoords); + std::vector * tgtCoords = & tgtCellCoords; + int tgtNbNodes = tgtCellCoords.size()/SPACEDIM; + NormalizedCellType tT=PlanarIntersector::_meshT.getTypeOfElement(OTT::indFC(icellT)); + bool isTargetQuad=CellModel::getCellModel(tT).isQuadratic(); + + typename MyMatrix::value_type& resRow=res[icellT]; + + // treat each source triangle + for(typename std::vector::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++) + { + int iS=*iter; + PlanarIntersector::getRealSourceCoordinates(OTT::indFC(iS),srcTriaCoords); + const ConnType *startOfCellNodeConn=PlanarIntersector::_connectS+OTT::conn2C(PlanarIntersector::_connIndexS[iS]); + if(SPACEDIM==3) + { + tgtCellCoordsTmp = tgtCellCoords; + tgtCoords = & tgtCellCoordsTmp; + orientation=PlanarIntersector::projectionThis(&tgtCellCoordsTmp[0], &srcTriaCoords[0], + tgtNbNodes, 3); + } + //double surf=orientation*intersectGeometryWithQuadrangle(quadrangle,targetCellCoordsTmp,isTargetQuad); + double surf=orientation*intersectGeoBary( *tgtCoords, isTargetQuad, &srcTriaCoords[0], nodeCeffs ); + surf=PlanarIntersector::getValueRegardingOption(surf); + if(surf!=0.) + { + for(int nodeIdS=0;nodeIdS<3;nodeIdS++) + { + ConnType curNodeS=startOfCellNodeConn[nodeIdS]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(curNodeS); + if(iterRes!=resRow.end()) + { + nodeCeffs[nodeIdS] += iterRes->second; + resRow.erase( curNodeS ); + } + resRow.insert(std::make_pair(curNodeS,nodeCeffs[nodeIdS])); + } + } + } + } +} +#endif diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.hxx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.hxx new file mode 100644 index 000000000..85c8c94fe --- /dev/null +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.hxx @@ -0,0 +1,63 @@ +// 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 __PolyhedronIntersectorP1P0Bary_HXX__ +#define __PolyhedronIntersectorP1P0Bary_HXX__ + +#include "Intersector3DP1P0Bary.hxx" +#include "SplitterTetra.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + + + /** + * \brief Class responsible for calculating intersection between a hexahedron target element and + * the source elements. + * + */ + template + class PolyhedronIntersectorP1P0Bary : public Intersector3DP1P0Bary + { + 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: + + PolyhedronIntersectorP1P0Bary(const MyMeshType& targetMesh, const MyMeshType& srcMesh, SplittingPolicy policy = GENERAL_24); + + ~PolyhedronIntersectorP1P0Bary(); + + void intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res); + + private: + void releaseArrays(); + private: + /// pointers to the SplitterTetra objects representing the tetrahedra + /// that result from the splitting of the hexahedron target cell + std::vector< SplitterTetra* > _tetra; + + SplitterTetra2 _split; + + }; +} + +#endif diff --git a/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx new file mode 100644 index 000000000..a44cb3b7f --- /dev/null +++ b/src/INTERP_KERNEL/PolyhedronIntersectorP1P0Bary.txx @@ -0,0 +1,141 @@ +// 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 __PolyhedronIntersectorP1P0Bary_TXX__ +#define __PolyhedronIntersectorP1P0Bary_TXX__ + +#include "PolyhedronIntersectorP1P0Bary.hxx" +#include "Intersector3DP1P0Bary.txx" +#include "MeshUtils.hxx" + +#include "SplitterTetra.txx" + +namespace INTERP_KERNEL +{ + + /** + * Constructor creating object from target cell global number + * The constructor first calculates the necessary nodes, + * (depending on the splitting policy) and then splits the hexahedron into + * tetrahedra, placing these in the internal vector _tetra. + * + * @param targetMesh mesh containing the target elements + * @param srcMesh mesh containing the source elements + * @param policy splitting policy to be used + * + * WARNING : in _split attribute, sourceMesh and targetMesh are switched in order to fit intersectCells feature. + */ + template + PolyhedronIntersectorP1P0Bary::PolyhedronIntersectorP1P0Bary(const MyMeshType& targetMesh, + const MyMeshType& srcMesh, + SplittingPolicy policy) + :Intersector3DP1P0Bary(targetMesh,srcMesh),_split(targetMesh,srcMesh,policy) + { + } + + /** + * Destructor. + * Liberates the SplitterTetra objects and potential sub-node points that have been allocated. + * + */ + template + PolyhedronIntersectorP1P0Bary::~PolyhedronIntersectorP1P0Bary() + { + releaseArrays(); + } + + template + void PolyhedronIntersectorP1P0Bary::releaseArrays() + { + for(typename std::vector< SplitterTetra* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + delete *iter; + _split.releaseArrays(); + _tetra.clear(); + } + + //================================================================================ + /*! + * \brief This method computes a value per each node of source cell for each target cell. + * \param srcCell - a source tetrahedron + * \param tgtCells - target elements + * \param res - matrix to fill in + */ + //================================================================================ + + template + void PolyhedronIntersectorP1P0Bary::intersectCells(ConnType tgtCell, + const std::vector& srcCells, + MyMatrix& res) + { + typename MyMatrix::value_type& resRow=res[tgtCell]; + + int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(tgtCell)); + releaseArrays(); + _split.splitTargetCell(tgtCell,nbOfNodesT,_tetra); + + for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) + { + // intersect a source tetrahedron with each target tetrahedron: get intersection volume and barycenter + double baryCentre[SPACEDIM], total_baryCentre[3] = { 0., 0., 0.}; + double interVolume = 0; + for(typename std::vector*>::iterator iterTetraT = _tetra.begin(); iterTetraT != _tetra.end(); ++iterTetraT) + { + SplitterTetra *tmp=*iterTetraT; + tmp->clearVolumesCache(); + double volume = tmp->intersectSourceCell(*iterCellS, baryCentre); + if ( volume > 0 ) + { + interVolume += volume; + for ( int i = 0; i < SPACEDIM; ++i ) + total_baryCentre[i] += baryCentre[i]*volume; + } + } + if(interVolume!=0) + { + for ( int i = 0; i < SPACEDIM; ++i ) + total_baryCentre[i] /= interVolume; + + // coordinates of the source tetrahedron + std::vector srcCellCoords(4); + for ( int n = 0; n < 4; ++n ) + srcCellCoords[ n ] = getCoordsOfNode( n, *iterCellS, Intersector3D::_src_mesh ); + + // compute barycentric coordinates + double baryCoords[4]; + barycentric_coords( srcCellCoords, total_baryCentre, baryCoords); + + // store coeffs of each node of the source tetrahedron + const ConnType *srcCellNodes=Intersector3D::_src_mesh.getConnectivityPtr()+OTT::conn2C(Intersector3D::_src_mesh.getConnectivityIndexPtr()[*iterCellS]); + for ( int n = 0; n < 4; ++n ) + { + double val = baryCoords[n] * interVolume; + ConnType curNodeS = srcCellNodes[n]; + typename MyMatrix::value_type::const_iterator iterRes=resRow.find(curNodeS); + if(iterRes!=resRow.end()) + { + val += iterRes->second; + resRow.erase( curNodeS ); + } + resRow.insert(std::make_pair(curNodeS,val)); + } + } + } + } +} + +#endif