From a61c29843f8bbc23d8437acf6248ca71e8fc3528 Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 16 Dec 2009 13:18:13 +0000 Subject: [PATCH] EHPOC Interpolators +InterpolationCU2D.hxx \ +InterpolationCU2D.txx \ --- src/INTERP_KERNEL/InterpolationCU2D.hxx | 50 ++++ src/INTERP_KERNEL/InterpolationCU2D.txx | 251 ++++++++++++++++++ src/INTERP_KERNEL/Makefile.am | 2 + .../Test/MEDCouplingBasicsTest.cxx | 79 +++++- .../Test/MEDCouplingBasicsTest.hxx | 3 + 5 files changed, 383 insertions(+), 2 deletions(-) create mode 100644 src/INTERP_KERNEL/InterpolationCU2D.hxx create mode 100644 src/INTERP_KERNEL/InterpolationCU2D.txx diff --git a/src/INTERP_KERNEL/InterpolationCU2D.hxx b/src/INTERP_KERNEL/InterpolationCU2D.hxx new file mode 100644 index 000000000..5ef80d4d4 --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationCU2D.hxx @@ -0,0 +1,50 @@ +// 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 +// +// File : InterpolationCU2D.hxx +// Created : Mon Dec 14 16:52:53 2009 +// Author : Edward AGAPOV (eap) + + +#ifndef __InterpolationCU2D_HXX__ +#define __InterpolationCU2D_HXX__ + +#include "Interpolation.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + class InterpolationCU2D : public Interpolation< InterpolationCU2D > + { + public: + InterpolationCU2D(); + InterpolationCU2D(const InterpolationOptions & io); + + template + int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method); + + template + int interpolateMeshes(const MyCMeshType& meshS, const MyUMeshType& meshT, MatrixType& result, const char *method); + + template + int interpolateMeshesRev(const MyUMeshType& meshS, const MyCMeshType& meshT, MatrixType& result, const char *method); + + }; +} + +#endif diff --git a/src/INTERP_KERNEL/InterpolationCU2D.txx b/src/INTERP_KERNEL/InterpolationCU2D.txx new file mode 100644 index 000000000..7dc991862 --- /dev/null +++ b/src/INTERP_KERNEL/InterpolationCU2D.txx @@ -0,0 +1,251 @@ +// 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 +// +// File : InterpolationCU2D.txx +// Created : Mon Dec 14 17:30:25 2009 +// Author : Edward AGAPOV (eap) + +#ifndef __InterpolationCU2D_TXX__ +#define __InterpolationCU2D_TXX__ + +#include "InterpolationCU2D.hxx" + +#include "Interpolation.txx" +#include "PlanarIntersector.txx" +#include "TriangulationIntersector.txx" + +#include + +// convert index "From Mesh Index" +#define _FMIU(i) OTT::ind2C((i)) +#define _FMIC(i) OTT::ind2C((i)) +// convert index "To Mesh Index" +#define _TMIU(i) OTT::indFC((i)) +#define _TMIC(i) OTT::indFC((i)) +// convert coord "From Mesh Coord" +#define _FMCOO(i) OTT::coo2C((i)) +// convert connectivity "From Mesh Connectivity" +#define _FMCON(i) OTT::conn2C((i)) + +namespace INTERP_KERNEL +{ + /** + * \defgroup InterpolationCU2D InterpolationCU2D + * \class InterpolationCU2D + * \brief Class used to calculate the volumes of intersection between the elements of a cartesian and an unstructured 2D meshes. + * + */ + //================================================================================ + /** + * Default constructor + * + */ + //================================================================================ + + InterpolationCU2D::InterpolationCU2D() + { + } + + InterpolationCU2D::InterpolationCU2D(const InterpolationOptions & io) + :Interpolation(io) + { + } + + //================================================================================ + /*! + * \brief A stub to anable inheriting from Interpolation + */ + //================================================================================ + + template + int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method) + { + throw INTERP_KERNEL::Exception("InterpolationCU2D class is intended to interpolate meshes of different nature: cartesian and unstructured"); + } + + //================================================================================ + /** + * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh. + * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the + * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the + * volume of intersection is calculated for the remaining pairs, and entered into the + * intersection matrix. + * + * The matrix is partially sparse : it is a vector of maps of integer - double pairs. + * It can also be an INTERP_KERNEL::Matrix object. + * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless + * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements + * which have a non-zero intersection volume with the target element. The vector has indices running from + * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however, + * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j]. + * + + * @param srcMesh 2-dimensional cartesian source mesh + * @param targetMesh 2-dimesional unstructured target mesh + * @param result matrix in which the result is stored + * @param method interpolation method + */ + //================================================================================ + + template + int InterpolationCU2D::interpolateMeshes(const MyCMeshType& src_mesh, + const MyUMeshType& tgt_mesh, + MatrixType& result, + const char * method) + { + typedef typename MyUMeshType::MyConnType ConnType; + + if ( std::string("P0P0") != method ) + throw Exception("Only P0P0 method is implemented so far"); + + // create empty maps for all target elements + result.resize( tgt_mesh.getNumberOfElements() ); + + const int ret = src_mesh.getNumberOfElements(); + + const double eps = getPrecision(); + const int dim = 2; + + const double* src_coords[ dim ]; + int src_nb_coords[ dim ]; + std::map< double, int> src_coord_to_index[ dim ]; + for ( int j = 0; j < dim; ++j ) + { + src_coords [j] = src_mesh.getCoordsAlongAxis( _TMIC( j )); + src_nb_coords[j] = src_mesh.nbCellsAlongAxis ( _TMIC( j )) + 1; + for (int i = 0; i < src_nb_coords[j]; ++i ) + src_coord_to_index[j].insert( std::make_pair( src_coords[j][i], i )); + } + + const unsigned long tgtu_nb_cells = tgt_mesh.getNumberOfElements(); + const ConnType * tgtu_conn = tgt_mesh.getConnectivityPtr(); + const ConnType * tgtu_conn_ind = tgt_mesh.getConnectivityIndexPtr(); + const double * tgtu_coord = tgt_mesh.getCoordinatesPtr(); + + double bb[2*dim]; + TriangulationIntersector + intersector (tgt_mesh,tgt_mesh, 0,0,0,0,0,0 ); + + // loop on unstructured tgt cells + + for(int iT=0; iT src_coords[0][ src_nb_coords[0]-1 ] - eps || + bb[2] > src_coords[1][ src_nb_coords[1]-1 ] - eps || + bb[1] < src_coords[0][0] + eps || + bb[3] < src_coords[1][0] + eps ) + continue; // no intersection + + // find structured src cells intersecting iT cell + int min_i[dim], max_i[dim]; + std::map< double, int>::iterator coo_ind; + for ( int j = 0; j < dim; ++j ) + { + coo_ind = src_coord_to_index[j].lower_bound( bb[2*j+1] - eps ); + if ( coo_ind == src_coord_to_index[j].end() ) + --coo_ind; + max_i[j] = coo_ind->second; + + coo_ind = src_coord_to_index[j].upper_bound( bb[2*j ] + eps ); + if ( coo_ind != src_coord_to_index[j].begin() ) + --coo_ind; + min_i[j] = coo_ind->second; + } + + // get tgt cell coords + std::vector tgtu_coords(2*tgtu_nb_nodes); + for (ConnType i=0; i eps ) + result[ iT ][ _TMIC( iX + iY * ( src_nb_coords[1]-1 )) ] = surf; + } + } + + return ret; + } + + //================================================================================ + /** + * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh. + * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the + * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the + * volume of intersection is calculated for the remaining pairs, and entered into the + * intersection matrix. + * + * The matrix is partially sparse : it is a vector of maps of integer - double pairs. + * It can also be an INTERP_KERNEL::Matrix object. + * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless + * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements + * which have a non-zero intersection volume with the target element. The vector has indices running from + * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however, + * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j]. + * + + * @param srcMesh 2-dimesional unstructured target mesh + * @param targetMesh 2-dimensional cartesian source mesh + * @param result matrix in which the result is stored + * @param method interpolation method + */ + //================================================================================ + + template + int InterpolationCU2D::interpolateMeshesRev(const MyUMeshType& meshS, const MyCMeshType& meshT, MatrixType& result, const char *method) + { + MatrixType revResult; + int sizeT = interpolateMeshes( meshT, meshS, revResult, method ); + int sizeS = revResult.size(); + result.resize( sizeT ); + + for ( int iS = 0; iS < sizeS; ++iS ) + { + typename MatrixType::value_type & row = revResult[iS]; + typename MatrixType::value_type::iterator iT_surf = row.begin(); + for ( ; iT_surf != row.end(); ++iT_surf ) + result[ iT_surf->first ][ iS ] = iT_surf->second; + } + return sizeS; + } + +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 6b6cfe576..77daecfa3 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -55,6 +55,8 @@ InterpolationPlanar.hxx \ InterpolationPlanar.txx \ InterpolationCC.hxx \ InterpolationCC.txx \ +InterpolationCU2D.hxx \ +InterpolationCU2D.txx \ InterpolationUtils.hxx \ Intersector3D.hxx \ Intersector3D.txx \ diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx index e84cea480..a42afbae0 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx @@ -24,6 +24,7 @@ #include "Interpolation3DSurf.txx" #include "Interpolation3D.txx" #include "InterpolationCC.txx" +#include #include "MEDCouplingNormalizedUnstructuredMesh.txx" #include "MEDCouplingNormalizedCartesianMesh.txx" @@ -2261,8 +2262,8 @@ void MEDCouplingBasicsTest::testInterpolationCC() mesh[i]->setCoords( coords, coords, coords ); coords->decrRef(); } - MEDCouplingNormalizedCartesianMesh<3,3> targetWrapper(mesh[1]); - MEDCouplingNormalizedCartesianMesh<3,3> sourceWrapper(mesh[0]); + MEDCouplingNormalizedCartesianMesh<3> targetWrapper(mesh[1]); + MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(mesh[0]); CPPUNIT_ASSERT_EQUAL( 27,int( sourceWrapper.getNumberOfElements())); CPPUNIT_ASSERT_EQUAL( 3, int( sourceWrapper.nbCellsAlongAxis(0))); CPPUNIT_ASSERT_EQUAL( 3, int( sourceWrapper.nbCellsAlongAxis(1))); @@ -2312,6 +2313,59 @@ void MEDCouplingBasicsTest::testInterpolationCC() mesh[1]->decrRef(); } +void MEDCouplingBasicsTest::testInterpolationCU2D() +{ + MEDCouplingCMesh* meshC = MEDCouplingCMesh::New(); + DataArrayDouble* coords = DataArrayDouble::New(); + double arr[4] = { 0/3, 1/3., 2/3., 3/3. }; + coords->useArray( arr, /*ownership=*/false, CPP_DEALLOC, 4, 1 ); + meshC->setCoords( coords, coords ); + coords->decrRef(); + + MEDCouplingUMesh * meshU = buildCU2DMesh_U(); + + MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(meshC); + MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(meshU); + INTERP_KERNEL::InterpolationCU2D myInterpolator; + vector > res; + myInterpolator.setPrecision(1e-12); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + +// cout.precision(18); +// for ( int i = 0; i < (int)res.size(); ++i ) +// for ( map::iterator s_v = res[i].begin(); s_v != res[i].end(); ++s_v) +// { +// cout << "CPPUNIT_ASSERT_DOUBLES_EQUAL( "<second<<" ,res["<first<<"],precis);"<decrRef(); + meshU->decrRef(); +} + void MEDCouplingBasicsTest::test2DInterpP0IntegralUniform() { MEDCouplingUMesh *targetMesh=build2DTargetMesh_1(); @@ -2791,6 +2845,27 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_2() return targetMesh; } +MEDCouplingUMesh *MEDCouplingBasicsTest::buildCU2DMesh_U() +{ + double coords[18]={0.0,0.0, 0.5,0.0, 1.0,0.0, 0.0,0.5, 0.5,0.5, 1.0,0.5, 0.0,1.0, 0.5,1.0, 1.0,1.0 }; + int conn[18]={0,1,4,3, 3,4,7,6, 4,5,8,7, 1,5,4, 1,2,5 }; + MEDCouplingUMesh *mesh=MEDCouplingUMesh::New(); + mesh->setMeshDimension(2); + mesh->allocateCells(5); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+4); + mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+8); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+12); + mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+15); + mesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,2); + std::copy(coords,coords+18,myCoords->getPointer()); + mesh->setCoords(myCoords); + myCoords->decrRef(); + return mesh; +} + MEDCouplingUMesh *MEDCouplingBasicsTest::build3DSurfSourceMesh_1() { double sourceCoords[12]={-0.3,-0.3,0.5, 0.7,-0.3,1.5, -0.3,0.7,0.5, 0.7,0.7,1.5}; diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 2bf2f3026..0f72805b5 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -73,6 +73,7 @@ namespace ParaMEDMEM CPPUNIT_TEST( test3DSurfInterpP0P0_2 ); CPPUNIT_TEST( test3DSurfInterpP0P0_3 ); CPPUNIT_TEST( testInterpolationCC ); + CPPUNIT_TEST( testInterpolationCU2D ); CPPUNIT_TEST( test3DInterpP0P0_1 ); CPPUNIT_TEST( test3DInterpP0P0PL_1 ); CPPUNIT_TEST( test3DInterpP0P0PL_2 ); @@ -149,6 +150,7 @@ namespace ParaMEDMEM void test3DInterpP1P1_1(); void test3DInterpP1P1PL_1(); void testInterpolationCC(); + void testInterpolationCU2D(); void test3DInterpP0P0Empty(); void test2DInterpP0IntegralUniform(); void test3DSurfInterpP0IntegralUniform(); @@ -167,6 +169,7 @@ namespace ParaMEDMEM MEDCouplingUMesh *build2DTargetMesh_1(); MEDCouplingUMesh *build2DTargetMeshPerm_1(); MEDCouplingUMesh *build2DTargetMesh_2(); + MEDCouplingUMesh *buildCU2DMesh_U(); MEDCouplingUMesh *build3DSurfSourceMesh_1(); MEDCouplingUMesh *build3DSurfSourceMesh_2(); MEDCouplingUMesh *build3DSurfTargetMesh_1(); -- 2.39.2