From: vbd Date: Fri, 1 Feb 2008 10:27:33 +0000 (+0000) Subject: *** empty log message *** X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=e880d9e0c8d3acd8574d4aad58b887052d712f66;p=tools%2Fmedcoupling.git *** empty log message *** --- diff --git a/src/INTERP_KERNEL/BBTree.H b/src/INTERP_KERNEL/BBTree.H index 211a15032..cc13b8166 100644 --- a/src/INTERP_KERNEL/BBTree.H +++ b/src/INTERP_KERNEL/BBTree.H @@ -28,7 +28,7 @@ private: public: -BBTree(double* bbs, int* elems, int level, int nbelems): +BBTree(double* bbs, int nbelems, int* elems=0, int level=0): _left(0), _right(0), _level(level), _bb(bbs), _terminal(false),_nbelems(nbelems) { if (nbelems < MIN_NB_ELEMS || level> MAX_LEVEL) @@ -89,8 +89,8 @@ BBTree(double* bbs, int* elems, int level, int nbelems): } _max_left=max_left; _min_right=min_right; - _left=new BBTree(bbs, &(new_elems_left[0]), level+1, new_elems_left.size()); - _right=new BBTree(bbs, &(new_elems_right[0]), level+1, new_elems_right.size()); + _left=new BBTree(bbs, new_elems_left.size(), &(new_elems_left[0]), level+1); + _right=new BBTree(bbs, new_elems_right.size(), &(new_elems_right[0]), level+1); } @@ -123,7 +123,7 @@ void getIntersectingElems(const double* bb, vector& elems) const bool intersects = true; for (int idim=0; idim1e-12|| bb_ptr[idim*2+1]-bb[idim*2]<-1e-12) + if (bb_ptr[idim*2]-bb[idim*2+1]>-1e-12|| bb_ptr[idim*2+1]-bb[idim*2]<1e-12) intersects=false; } if (intersects) diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx index a5fdc9869..39a4e6063 100644 --- a/src/INTERP_KERNEL/Interpolation.hxx +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -6,7 +6,7 @@ #include #include "MEDMEM_Mesh.hxx" - +#include "MEDMEM_Matrix.hxx" /** * \mainpage * Status : documentation of 3D - part of intersection matrix calculation more or less complete @@ -24,9 +24,7 @@ namespace MEDMEM Interpolation(){} virtual ~Interpolation(){} - //interpolation of two triangular meshes. - virtual std::vector > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0; - + virtual void printInfo(int level)=0; }; }; diff --git a/src/INTERP_KERNEL/Interpolation2D.cxx b/src/INTERP_KERNEL/Interpolation2D.cxx index 99861350e..a4f5336b8 100755 --- a/src/INTERP_KERNEL/Interpolation2D.cxx +++ b/src/INTERP_KERNEL/Interpolation2D.cxx @@ -1,14 +1,7 @@ #include "MEDMEM_Mesh.hxx" #include "Interpolation2D.hxx" -#include "PlanarIntersector.hxx" -#include "PlanarIntersector.H" -#include "TriangulationIntersector.hxx" -#include "TriangulationIntersector.cxx" -#include "ConvexIntersector.hxx" -#include "ConvexIntersector.cxx" -#include "GenericIntersector.hxx" -#include "GenericIntersector.cxx" -#include "BBTree.H" +#include "Interpolation2D.txx" + #include using namespace MED_EN; @@ -53,183 +46,358 @@ namespace MEDMEM _PrintLevel=PrintLevel; _Intersection_type=intersection_type; } +// std::vector > Interpolation2D::interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2) +// { +// std::vector > matrix; +// interpolateMeshes(mesh1,mesh2,matrix); +// return matrix; +// } +// +// MEDMEM::MEDMEM_Matrix Interpolation2D::interpolate(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2) +// { +// MEDMEM::MEDMEM_Matrix matrix; +// interpolateMeshes(mesh1,mesh2,matrix); +// return matrix; +// } - - /** \brief Main function to interpolate triangular or quadrangular meshes. - \details The algorithm proceeds 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 by an object of type Intersector2D 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. - * 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 (#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 myMesh_S 2D source mesh - * @param myMesh_P 2D target mesh - * @return vector containing for each element i of the source mesh, a map giving for each element j - * of the target mesh which i intersects, the area of the intersection - * - */ - - vector > Interpolation2D::interpolateMeshes(const MEDMEM::MESH& myMesh_S, - const MEDMEM::MESH& myMesh_P) - { - long global_start =clock(); - int counter=0; - /***********************************************************/ - /* Check both meshes are made of triangles and quadrangles */ - /***********************************************************/ - - int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL); - int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL); - if ( NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2) - cout<<"Interpolation2D::intersectMeshes: one of the two meshes contains more than two types of cells"< > BoxS=myMesh_S.getBoundingBox(); - vector > BoxP=myMesh_P.getBoundingBox(); - double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) - +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])); - double DimCaracteristic_S=diagonalS/nbMaille_S; - double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) - +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])); - double DimCaracteristic_P=diagonalP/nbMaille_P; - - _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); - if (_PrintLevel>=1) - { - cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, - 0, _PrintLevel); - break; - case MEDMEM::Convex: - intersector=new ConvexIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, - 0, 0, _PrintLevel); - break; - case MEDMEM::Generic: - intersector=new GenericIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, - 0, 0, _PrintLevel); - break; - } - - /****************************************************************/ - /* Create a search tree based on the bounding boxes */ - /* Instanciate the intersector and initialise the result vector */ - /****************************************************************/ - - long start_filtering=clock(); - - vector bbox; - intersector->createBoundingBoxes<2>(myMesh_S,bbox); // create the bounding boxes - BBTree<2> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure +// +// +// /** \brief Main function to interpolate triangular or quadrangular meshes. +// \details The algorithm proceeds 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 by an object of type Intersector2D 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. +// * 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 (#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 myMesh_S 2D source mesh +// * @param myMesh_P 2D target mesh +// * @return vector containing for each element i of the source mesh, a map giving for each element j +// * of the target mesh which i intersects, the area of the intersection +// * +// */ +//template +// void Interpolation2D::interpolateMeshes(const MEDMEM::MESH& myMesh_S, +// const MEDMEM::MESH& myMesh_P, +// T& matrix) +// { +// long global_start =clock(); +// int counter=0; +// /***********************************************************/ +// /* Check both meshes are made of triangles and quadrangles */ +// /***********************************************************/ +// +// int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL); +// int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL); +// if ( NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2) +// cout<<"Interpolation2D::intersectMeshes: one of the two meshes contains more than two types of cells"< > BoxS=myMesh_S.getBoundingBox(); +// vector > BoxP=myMesh_P.getBoundingBox(); +// double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) +// +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])); +// double DimCaracteristic_S=diagonalS/nbMaille_S; +// double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) +// +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])); +// double DimCaracteristic_P=diagonalP/nbMaille_P; +// +// _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); +// if (_PrintLevel>=1) +// { +// cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// 0, _PrintLevel); +// break; +// case MEDMEM::Convex: +// intersector=new ConvexIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, +// 0, 0, _PrintLevel); +// break; +// case MEDMEM::Generic: +// intersector=new GenericIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// 0, 0, _PrintLevel); +// break; +// } +// +// /****************************************************************/ +// /* Create a search tree based on the bounding boxes */ +// /* Instanciate the intersector and initialise the result vector */ +// /****************************************************************/ +// +// long start_filtering=clock(); +// +// vector bbox; +// intersector->createBoundingBoxes<2>(myMesh_S,bbox); // create the bounding boxes +// BBTree<2> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure +// +// long end_filtering=clock(); +// +// map MAP_init; +// vector > result_areas(nbMaille_P,MAP_init);//on initialise. +// +// /****************************************************/ +// /* Loop on the target cells - core of the algorithm */ +// /****************************************************/ +// const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL); +// int type_max_index=0;//maximum cell number for a given type +// int type_min_index=0;//minimum cell number for a given type +// int i_P=0;//global index of cell +// +// long start_intersection=clock(); +// +// for (int itype = 0; itype intersecting_elems; +// double bb[4]; +// intersector->getElemBB<2>(bb,myMesh_P,i_P+1,nb_nodesP); +// my_tree.getIntersectingElems(bb, intersecting_elems); +// int nb_intersecting_elems = intersecting_elems.size(); +// //browsing all the i_S (from mesh S) elems that can +// //intersect elem i_P (from mesh P) +// for (int ielem=0; ielemintersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); +// (result_areas[i_P]).insert(make_pair(i_S+1,surf)); +// counter++; +// //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); +// } +// intersecting_elems.clear(); +// } +// type_min_index = type_max_index; +// } +// delete intersector; +// +// /***********************************/ +// /* DEBUG prints */ +// /***********************************/ +// +// if (_PrintLevel >=1) +// { +// long end_intersection=clock(); +// if (_PrintLevel >=2) +// { +// cout<::iterator surface; +// total_interm=0.0; +// for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++) +// { +// cout<<" ("< > BoxS=myMesh_S.getBoundingBox(); +// vector > BoxP=myMesh_P.getBoundingBox(); +// double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) +// +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])); +// double DimCaracteristic_S=diagonalS/nbMaille_S; +// double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) +// +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])); +// double DimCaracteristic_P=diagonalP/nbMaille_P; +// +// _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); +// if (_PrintLevel>=1) +// { +// cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// 0, _PrintLevel); +// break; +// case MEDMEM::Convex: +// intersector=new ConvexIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, +// 0, 0, _PrintLevel); +// break; +// case MEDMEM::Generic: +// intersector=new GenericIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// 0, 0, _PrintLevel); +// break; +// } +// +// /****************************************************************/ +// /* Create a search tree based on the bounding boxes */ +// /* Instanciate the intersector and initialise the result vector */ +// /****************************************************************/ +// +// long start_filtering=clock(); +// +// vector bbox; +// intersector->createBoundingBoxes<2>(myMesh_S,bbox); // create the bounding boxes +// BBTree<2> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure +// +// long end_filtering=clock(); +// +// map MAP_init; +// MEDMEM_Matrix* matrix=new MEDMEM_Matrix(nbMaille_P); +// +// /****************************************************/ +// /* Loop on the target cells - core of the algorithm */ +// /****************************************************/ +// const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL); +// int type_max_index=0;//maximum cell number for a given type +// int type_min_index=0;//minimum cell number for a given type +// int i_P=0;//global index of cell +// +// long start_intersection=clock(); +// +// for (int itype = 0; itype intersecting_elems; +// double bb[4]; +// intersector->getElemBB<2>(bb,myMesh_P,i_P+1,nb_nodesP); +// my_tree.getIntersectingElems(bb, intersecting_elems); +// int nb_intersecting_elems = intersecting_elems.size(); +// //browsing all the i_S (from mesh S) elems that can +// //intersect elem i_P (from mesh P) +// for (int ielem=0; ielemintersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); +// matrix->setIJ(i_P+1,i_S+1,surf); +// //(result_areas[i_P]).insert(make_pair(i_S+1,surf)); +// counter++; +// //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); +// } +// intersecting_elems.clear(); +// } +// type_min_index = type_max_index; +// } +// delete intersector; +// +// /***********************************/ +// /* DEBUG prints */ +// /***********************************/ +// +// if (_PrintLevel >=1) +// { +// long end_intersection=clock(); +// if (_PrintLevel >=2) +// { +// cout<<"Printing intersection matrix"; +// cout<<*matrix; +// } +// +// cout<< "Filtering time= " << end_filtering-start_filtering<< endl; +// cout<< "Intersection time= " << end_intersection-start_intersection<< endl; +// long global_end =clock(); +// cout<< "Number of computed intersections = " << counter << endl; +// cout<< "Global time= " << global_end - global_start << endl; +// } +// +// return matrix; +// } - long end_filtering=clock(); - - map MAP_init; - vector > result_areas(nbMaille_P,MAP_init);//on initialise. - - /****************************************************/ - /* Loop on the target cells - core of the algorithm */ - /****************************************************/ - const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL); - int type_max_index=0;//maximum cell number for a given type - int type_min_index=0;//minimum cell number for a given type - int i_P=0;//global index of cell - - long start_intersection=clock(); - - for (int itype = 0; itype intersecting_elems; - double bb[4]; - intersector->getElemBB<2>(bb,myMesh_P,i_P+1,nb_nodesP); - my_tree.getIntersectingElems(bb, intersecting_elems); - int nb_intersecting_elems = intersecting_elems.size(); - //browsing all the i_S (from mesh S) elems that can - //intersect elem i_P (from mesh P) - for (int ielem=0; ielemintersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); - (result_areas[i_P]).insert(make_pair(i_S+1,surf)); - counter++; - //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); - } - intersecting_elems.clear(); - } - type_min_index = type_max_index; - } - delete intersector; - - /***********************************/ - /* DEBUG prints */ - /***********************************/ - - if (_PrintLevel >=1) - { - long end_intersection=clock(); - if (_PrintLevel >=2) - { - cout<::iterator surface; - total_interm=0.0; - for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++) - { - cout<<" ("< + void Interpolation2D::interpolateMeshes(const MEDMEM::MESH& myMesh_S, + const MEDMEM::MESH& myMesh_P, + T& matrix) + { + long global_start =clock(); + int counter=0; + /***********************************************************/ + /* Check both meshes are made of triangles and quadrangles */ + /***********************************************************/ + + int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL); + int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL); + if ( NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2) + cout<<"Interpolation2D::intersectMeshes: one of the two meshes contains more than two types of cells"< > BoxS=myMesh_S.getBoundingBox(); + vector > BoxP=myMesh_P.getBoundingBox(); + double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) + +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1])); + double DimCaracteristic_S=diagonalS/nbMaille_S; + double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) + +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1])); + double DimCaracteristic_P=diagonalP/nbMaille_P; + + _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); + if (_PrintLevel>=1) + { + cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, + 0, _PrintLevel); + break; + case MEDMEM::Convex: + intersector=new ConvexIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, + 0, 0, _PrintLevel); + break; + case MEDMEM::Generic: + intersector=new GenericIntersector<2>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, + 0, 0, _PrintLevel); + break; + } + + /****************************************************************/ + /* Create a search tree based on the bounding boxes */ + /* Instanciate the intersector and initialise the result vector */ + /****************************************************************/ + + long start_filtering=clock(); + + vector bbox; + intersector->createBoundingBoxes<2>(myMesh_S,bbox); // create the bounding boxes + BBTree<2> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure + + long end_filtering=clock(); + + matrix.resize(nbMaille_P); + + /****************************************************/ + /* Loop on the target cells - core of the algorithm */ + /****************************************************/ + const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL); + int type_max_index=0;//maximum cell number for a given type + int type_min_index=0;//minimum cell number for a given type + int i_P=0;//global index of cell + + long start_intersection=clock(); + + for (int itype = 0; itype intersecting_elems; + double bb[4]; + intersector->getElemBB<2>(bb,myMesh_P,i_P+1,nb_nodesP); + my_tree.getIntersectingElems(bb, intersecting_elems); + int nb_intersecting_elems = intersecting_elems.size(); + //browsing all the i_S (from mesh S) elems that can + //intersect elem i_P (from mesh P) + for (int ielem=0; ielemintersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); + (matrix[i_P]).insert(make_pair(i_S+1,surf)); + counter++; + //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); + } + intersecting_elems.clear(); + } + type_min_index = type_max_index; + } + delete intersector; + + /***********************************/ + /* DEBUG prints */ + /***********************************/ + + if (_PrintLevel >=1) + { + long end_intersection=clock(); +// if (_PrintLevel >=2) +// { +// cout<::iterator surface; +// total_interm=0.0; +// for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++) +// { +// cout<<" ("< #else // use BBTree class - #include "BBTree.H" #endif -namespace MEDMEM -{ - /** - * \defgroup interpolation3D Interpolation3D - * \class Interpolation3D - * \brief Class used to calculate the volumes of intersection between the elements of two 3D meshes. - * - */ - /** - * Default constructor - * - */ - Interpolation3D::Interpolation3D() - { - } - - /** - * Destructor - * - */ - Interpolation3D::~Interpolation3D() - { - } - - /** - * Main interface method of the class, which verifies that the meshes are valid for the calculation, - * and then returns the matrix of intersection volumes. - * - * @param srcMesh 3-dimensional source mesh - * @param targetMesh 3-dimesional target mesh, containing only tetraedra - * @return vector containing for each element i of the source mesh, a map giving for each element j - * of the target mesh which i intersects, the volume of the intersection - */ - IntersectionMatrix Interpolation3D::interpolateMeshes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) - { - // it seems wasteful to make a copy here - IntersectionMatrix matrix; - - // we should do more sanity checking here - eliminating meshes that - // cannot be treated +namespace MEDMEM { +/** + * \defgroup interpolation3D Interpolation3D + * \class Interpolation3D + * \brief Class used to calculate the volumes of intersection between the elements of two 3D meshes. + * + */ +/** + * Default constructor + * + */ +Interpolation3D::Interpolation3D() { +} + +/** + * Destructor + * */ +Interpolation3D::~Interpolation3D() { +} + +/** + * Main interface method of the class, which verifies that the meshes are valid for the calculation, + * and then returns the matrix of intersection volumes. + * + * @param srcMesh 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @return vector containing for each element i of the source mesh, a map giving for each element j + * of the target mesh which i intersects, the volume of the intersection + */ +MEDMEM::MEDMEM_Matrix* Interpolation3D::interpolate( + const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) { + // it seems wasteful to make a copy here + //IntersectionMatrix matrix + MEDMEM::MEDMEM_Matrix* matrix; + // we should do more sanity checking here - eliminating meshes that + // cannot be treated + + calculateIntersectionVolumes(srcMesh, targetMesh, matrix); + return matrix; +} + +IntersectionMatrix Interpolation3D::interpolateMeshes( + const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) { + // it seems wasteful to make a copy here + IntersectionMatrix matrix; - calculateIntersectionVolumes(srcMesh, targetMesh, matrix); - return matrix; - } - - /** - * 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 by an object of type Intersector3D 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. - * 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 (#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 3-dimensional source mesh - * @param targetMesh 3-dimesional target mesh, containing only tetraedra - * @param matrix vector of maps in which the result is stored - * - */ - void Interpolation3D::calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix) { - - // create MeshElement objects corresponding to each element of the two meshes - const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); - const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS); - - LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); - - MeshElement* srcElems[numSrcElems]; - MeshElement* targetElems[numTargetElems]; - - std::map indices; - - for(int i = 0 ; i < numSrcElems ; ++i) - { - //const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); - srcElems[i] = new MeshElement(i + 1, srcMesh); - } - - for(int i = 0 ; i < numTargetElems ; ++i) - { - //const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); - targetElems[i] = new MeshElement(i + 1, targetMesh); - } - - // create empty maps for all source elements - matrix.resize(numTargetElems); + // we should do more sanity checking here - eliminating meshes that + // cannot be treated + + calculateIntersectionVolumes(srcMesh, targetMesh, matrix); + return matrix; +} + +/** + * 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 by an object of type Intersector3D 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. + * 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 (#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 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param matrix vector of maps in which the result is stored + * + */ +void Interpolation3D::calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, + const MEDMEM::MESH& targetMesh, MEDMEM::MEDMEM_Matrix*& matrix) { + + // create MeshElement objects corresponding to each element of the two meshes + const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + MeshElement* srcElems[numSrcElems]; + MeshElement* targetElems[numTargetElems]; + + std::map indices; + + for (int i = 0; i < numSrcElems; ++i) { + //const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); + srcElems[i] = new MeshElement(i + 1, srcMesh); + } + + for (int i = 0; i < numTargetElems; ++i) { + //const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); + targetElems[i] = new MeshElement(i + 1, targetMesh); + } + + // create empty maps for all source elements + matrix=new MEDMEM::MEDMEM_Matrix(numTargetElems); #ifdef USE_RECURSIVE_BBOX_FILTER - - /* - * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh - * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming - * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the - * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a - * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements - * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh. - */ - - // create initial RegionNode and fill up its source region with all the source mesh elements and - // its target region with all the target mesh elements whose bounding box - // intersects that of the source region - - RegionNode* firstNode = new RegionNode(); - - MeshRegion& srcRegion = firstNode->getSrcRegion(); - - for(int i = 0 ; i < numSrcElems ; ++i) - { - srcRegion.addElement(srcElems[i], srcMesh); - } - - MeshRegion& targetRegion = firstNode->getTargetRegion(); - - for(int i = 0 ; i < numTargetElems ; ++i) - { - if(!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]) )) - { - targetRegion.addElement(targetElems[i], targetMesh); - } - } - - // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and - // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the - // elements of the target mesh whose bounding box intersects the corresponding part - // Continue until the source region contains only one element, at which point the intersection volumes are - // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. - - stack nodes; - nodes.push(firstNode); - - while(!nodes.empty()) - { - RegionNode* currNode = nodes.top(); - nodes.pop(); - LOG(4, "Popping node "); - - if(currNode->getTargetRegion().getNumberOfElements() == 1) - { - // calculate volumes - LOG(4, " - One element"); - - MeshElement* targetElement = *(currNode->getTargetRegion().getBeginElements()); - - // NB : srcElement indices are from 0 .. numSrcElements - 1 - // targetElement indicies from 1 .. numTargetElements - // maybe this is not ideal ... - const int targetIdx = targetElement->getIndex(); - - - TargetIntersector* intersector; - - switch(targetElement->getType()) - { - case MED_TETRA4: - intersector = new IntersectorTetra(srcMesh, targetMesh, targetIdx); - break; - - case MED_HEXA8: - intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx); - break; - - default: - assert(false); - } - - - - for(vector::const_iterator iter = currNode->getSrcRegion().getBeginElements() ; - iter != currNode->getSrcRegion().getEndElements() ; ++iter) - { - - const int srcIdx = (*iter)->getIndex(); - const double vol = intersector->intersectSourceCell(srcIdx); - - if(vol != 0.0) - { - matrix[targetIdx - 1][srcIdx] = vol; - LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); - - } - } - - delete intersector; - - } - else // recursion - { - - LOG(4, " - Recursion"); - - RegionNode* leftNode = new RegionNode(); - RegionNode* rightNode = new RegionNode(); - - // split current source region - //} decide on axis - static BoundingBox::BoxCoord axis = BoundingBox::XMAX; - - currNode->getTargetRegion().split(leftNode->getTargetRegion(), rightNode->getTargetRegion(), axis, targetMesh); - - LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements() - << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements() - << " elements"); - - // ugly hack to avoid problem with enum which does not start at 0 - // I guess I ought to implement ++ for it instead ... - // Anyway, it basically chooses the next axis, cyclically - axis = (axis != BoundingBox::ZMAX) ? static_cast(axis + 1) : BoundingBox::XMAX; - - // add source elements of current node that overlap the target regions of the new nodes - LOG(5, " -- Adding source elements"); - int numLeftElements = 0; - int numRightElements = 0; - for(vector::const_iterator iter = currNode->getSrcRegion().getBeginElements() ; - iter != currNode->getSrcRegion().getEndElements() ; ++iter) - { - LOG(6, " --- New target node"); - - if(!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) - { - leftNode->getSrcRegion().addElement(*iter, srcMesh); - ++numLeftElements; - } - - if(!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) - { - rightNode->getSrcRegion().addElement(*iter, srcMesh); - ++numRightElements; - } - - } - - LOG(5, "Left src region has " << numLeftElements << " elements and right src region has " - << numRightElements << " elements"); - - // push new nodes on stack - if(numLeftElements != 0) - { - nodes.push(leftNode); - } - else - { - delete leftNode; - } - - if(numRightElements != 0) - { - nodes.push(rightNode); - } - else - { - delete rightNode; - } - } - - // all nodes are deleted here - delete currNode; - - LOG(4, "Next iteration. Nodes left : " << nodes.size()); - } - - + + /* + * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh + * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming + * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the + * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a + * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements + * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh. + */ + + // create initial RegionNode and fill up its source region with all the source mesh elements and + // its target region with all the target mesh elements whose bounding box + // intersects that of the source region + + RegionNode* firstNode = new RegionNode(); + + MeshRegion& srcRegion = firstNode->getSrcRegion(); + + for (int i = 0; i < numSrcElems; ++i) { + srcRegion.addElement(srcElems[i], srcMesh); + } + + MeshRegion& targetRegion = firstNode->getTargetRegion(); + + for (int i = 0; i < numTargetElems; ++i) { + if (!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]))) { + targetRegion.addElement(targetElems[i], targetMesh); + } + } + + // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and + // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the + // elements of the target mesh whose bounding box intersects the corresponding part + // Continue until the source region contains only one element, at which point the intersection volumes are + // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. + + stack nodes; + nodes.push(firstNode); + + while (!nodes.empty()) { + RegionNode* currNode = nodes.top(); + nodes.pop();LOG(4, "Popping node "); + + if (currNode->getTargetRegion().getNumberOfElements() == 1) { + // calculate volumes + LOG(4, " - One element"); + + MeshElement* targetElement = *(currNode->getTargetRegion().getBeginElements()); + + // NB : srcElement indices are from 0 .. numSrcElements - 1 + // targetElement indicies from 1 .. numTargetElements + // maybe this is not ideal ... + const int targetIdx = targetElement->getIndex(); + + TargetIntersector* intersector; + + switch (targetElement->getType()) { + case MED_TETRA4: + intersector = new IntersectorTetra(srcMesh, targetMesh, targetIdx); + break; + + case MED_HEXA8: + intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx); + break; + + default: + assert(false); + } + + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { + + const int srcIdx = (*iter)->getIndex(); + const double vol = intersector->intersectSourceCell(srcIdx); + + if (vol != 0.0) { + matrix->setIJ(targetIdx,srcIdx,vol); + //matrix[targetIdx - 1][srcIdx] = vol; + LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix->getIJ(srcIdx,targetIdx)); + + //LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); + + } + } + + delete intersector; + + } else // recursion + { + + LOG(4, " - Recursion"); + + RegionNode* leftNode = new RegionNode(); + RegionNode* rightNode = new RegionNode(); + + // split current source region + //} decide on axis + static BoundingBox::BoxCoord axis = BoundingBox::XMAX; + + currNode->getTargetRegion().split(leftNode->getTargetRegion(), + rightNode->getTargetRegion(), axis, targetMesh); + + LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements() + << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements() + << " elements"); + + // ugly hack to avoid problem with enum which does not start at 0 + // I guess I ought to implement ++ for it instead ... + // Anyway, it basically chooses the next axis, cyclically + axis + = (axis != BoundingBox::ZMAX) ? static_cast(axis + + 1) + : BoundingBox::XMAX; + + // add source elements of current node that overlap the target regions of the new nodes + LOG(5, " -- Adding source elements"); + int numLeftElements = 0; + int numRightElements = 0; + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { + LOG(6, " --- New target node"); + + if (!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + leftNode->getSrcRegion().addElement(*iter, srcMesh); + ++numLeftElements; + } + + if (!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + rightNode->getSrcRegion().addElement(*iter, srcMesh); + ++numRightElements; + } + + } + + LOG(5, "Left src region has " << numLeftElements << " elements and right src region has " + << numRightElements << " elements"); + + // push new nodes on stack + if (numLeftElements != 0) { + nodes.push(leftNode); + } else { + delete leftNode; + } + + if (numRightElements != 0) { + nodes.push(rightNode); + } else { + delete rightNode; + } + } + + // all nodes are deleted here + delete currNode; + + LOG(4, "Next iteration. Nodes left : " << nodes.size()); + } #else // Use BBTree - - // create BBTree structure - // - get bounding boxes - double bboxes[6 * numSrcElems]; - int srcElemIdx[numSrcElems]; - for(int i = 0; i < numSrcElems ; ++i) - { - // get source bboxes in right order - const BoundingBox* box = srcElems[i]->getBoundingBox(); - bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); - bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); - bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); - bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); - bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); - bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); - - // source indices have to begin with zero for BBox, I think - srcElemIdx[i] = srcElems[i]->getIndex() - 1; - } - - BBTree<3> tree(bboxes, srcElemIdx, 0, numSrcElems); - - // for each target element, get source elements with which to calculate intersection - // - calculate intersection by calling intersectCells - for(int i = 0; i < numTargetElems; ++i) - { - const BoundingBox* box = targetElems[i]->getBoundingBox(); - const int targetIdx = targetElems[i]->getIndex(); - - // get target bbox in right order - double targetBox[6]; - targetBox[0] = box->getCoordinate(BoundingBox::XMIN); - targetBox[1] = box->getCoordinate(BoundingBox::XMAX); - targetBox[2] = box->getCoordinate(BoundingBox::YMIN); - targetBox[3] = box->getCoordinate(BoundingBox::YMAX); - targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); - targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); - - vector intersectElems; - - tree.getIntersectingElems(targetBox, intersectElems); - - // create intersector - IntersectorTetra intersector(srcMesh, targetMesh, targetIdx); - - for(vector::const_iterator iter = intersectElems.begin() ; iter != intersectElems.end() ; ++iter) - { - - const int srcIdx = *iter + 1; - const double vol = intersector.intersectSourceCell(srcIdx); - - if(vol != 0.0) - { - matrix[targetIdx - 1].insert(make_pair(srcIdx, vol)); - } - - } - } - + // create BBTree structure + // - get bounding boxes + double bboxes[6 * numSrcElems]; + int srcElemIdx[numSrcElems]; + for(int i = 0; i < numSrcElems; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + + // source indices have to begin with zero for BBox, I think + srcElemIdx[i] = srcElems[i]->getIndex() - 1; + } + + BBTree<3> tree(bboxes, srcElemIdx, 0, numSrcElems); + + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(int i = 0; i < numTargetElems; ++i) + { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + + vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + // create intersector + IntersectorTetra intersector(srcMesh, targetMesh, targetIdx); + + for(vector::const_iterator iter = intersectElems.begin(); iter != intersectElems.end(); ++iter) + { + + const int srcIdx = *iter + 1; + const double vol = intersector.intersectSourceCell(srcIdx); + + if(vol != 0.0) + { + matrix->setIJ(targetIdx,srcIdx,vol); + //matrix[targetIdx - 1].insert(make_pair(srcIdx, vol)); + } + + } + } + #endif + // free allocated memory + for (int i = 0; i < numSrcElems; ++i) { + delete srcElems[i]; + } + for (int i = 0; i < numTargetElems; ++i) { + delete targetElems[i]; + } + +} + +/** + * 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 by an object of type Intersector3D 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. + * 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 (#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 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param matrix vector of maps in which the result is stored + * + */ +void Interpolation3D::calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, + const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix) { + + // create MeshElement objects corresponding to each element of the two meshes + const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + MeshElement* srcElems[numSrcElems]; + MeshElement* targetElems[numTargetElems]; + + std::map indices; + + for (int i = 0; i < numSrcElems; ++i) { + //const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); + srcElems[i] = new MeshElement(i + 1, srcMesh); + } + + for (int i = 0; i < numTargetElems; ++i) { + //const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); + targetElems[i] = new MeshElement(i + 1, targetMesh); + } + + // create empty maps for all source elements + matrix.resize(numTargetElems); + +#ifdef USE_RECURSIVE_BBOX_FILTER + + /* + * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh + * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming + * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the + * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a + * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements + * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh. + */ + + // create initial RegionNode and fill up its source region with all the source mesh elements and + // its target region with all the target mesh elements whose bounding box + // intersects that of the source region + + RegionNode* firstNode = new RegionNode(); + + MeshRegion& srcRegion = firstNode->getSrcRegion(); + + for (int i = 0; i < numSrcElems; ++i) { + srcRegion.addElement(srcElems[i], srcMesh); + } + + MeshRegion& targetRegion = firstNode->getTargetRegion(); + + for (int i = 0; i < numTargetElems; ++i) { + if (!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]))) { + targetRegion.addElement(targetElems[i], targetMesh); + } + } + + // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and + // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the + // elements of the target mesh whose bounding box intersects the corresponding part + // Continue until the source region contains only one element, at which point the intersection volumes are + // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. + + stack nodes; + nodes.push(firstNode); + + while (!nodes.empty()) { + RegionNode* currNode = nodes.top(); + nodes.pop();LOG(4, "Popping node "); + + if (currNode->getTargetRegion().getNumberOfElements() == 1) { + // calculate volumes + LOG(4, " - One element"); + + MeshElement* targetElement = *(currNode->getTargetRegion().getBeginElements()); + + // NB : srcElement indices are from 0 .. numSrcElements - 1 + // targetElement indicies from 1 .. numTargetElements + // maybe this is not ideal ... + const int targetIdx = targetElement->getIndex(); + + TargetIntersector* intersector; + + switch (targetElement->getType()) { + case MED_TETRA4: + intersector = new IntersectorTetra(srcMesh, targetMesh, targetIdx); + break; + + case MED_HEXA8: + intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx); + break; + + default: + assert(false); + } + + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { - // free allocated memory - for(int i = 0 ; i < numSrcElems ; ++i) - { - delete srcElems[i]; - } - for(int i = 0 ; i < numTargetElems ; ++i) - { - delete targetElems[i]; - } + const int srcIdx = (*iter)->getIndex(); + const double vol = intersector->intersectSourceCell(srcIdx); + if (vol != 0.0) { + matrix[targetIdx - 1][srcIdx] = vol; + LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); - } + } + } + + delete intersector; + + } else // recursion + { + + LOG(4, " - Recursion"); + + RegionNode* leftNode = new RegionNode(); + RegionNode* rightNode = new RegionNode(); + + // split current source region + //} decide on axis + static BoundingBox::BoxCoord axis = BoundingBox::XMAX; + + currNode->getTargetRegion().split(leftNode->getTargetRegion(), + rightNode->getTargetRegion(), axis, targetMesh); + + LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements() + << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements() + << " elements"); + + // ugly hack to avoid problem with enum which does not start at 0 + // I guess I ought to implement ++ for it instead ... + // Anyway, it basically chooses the next axis, cyclically + axis + = (axis != BoundingBox::ZMAX) ? static_cast(axis + + 1) + : BoundingBox::XMAX; + + // add source elements of current node that overlap the target regions of the new nodes + LOG(5, " -- Adding source elements"); + int numLeftElements = 0; + int numRightElements = 0; + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { + LOG(6, " --- New target node"); + + if (!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + leftNode->getSrcRegion().addElement(*iter, srcMesh); + ++numLeftElements; + } + + if (!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + rightNode->getSrcRegion().addElement(*iter, srcMesh); + ++numRightElements; + } + + } + + LOG(5, "Left src region has " << numLeftElements << " elements and right src region has " + << numRightElements << " elements"); + + // push new nodes on stack + if (numLeftElements != 0) { + nodes.push(leftNode); + } else { + delete leftNode; + } + + if (numRightElements != 0) { + nodes.push(rightNode); + } else { + delete rightNode; + } + } + + // all nodes are deleted here + delete currNode; + + LOG(4, "Next iteration. Nodes left : " << nodes.size()); + } + +#else // Use BBTree + // create BBTree structure + // - get bounding boxes + double bboxes[6 * numSrcElems]; + int srcElemIdx[numSrcElems]; + for(int i = 0; i < numSrcElems; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + + // source indices have to begin with zero for BBox, I think + srcElemIdx[i] = srcElems[i]->getIndex() - 1; + } + + BBTree<3> tree(bboxes, srcElemIdx, 0, numSrcElems); + + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(int i = 0; i < numTargetElems; ++i) + { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + + vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + // create intersector + IntersectorTetra intersector(srcMesh, targetMesh, targetIdx); + + for(vector::const_iterator iter = intersectElems.begin(); iter != intersectElems.end(); ++iter) + { + + const int srcIdx = *iter + 1; + const double vol = intersector.intersectSourceCell(srcIdx); + + if(vol != 0.0) + { + matrix[targetIdx - 1].insert(make_pair(srcIdx, vol)); + } + + } + } + +#endif + + // free allocated memory + for (int i = 0; i < numSrcElems; ++i) { + delete srcElems[i]; + } + for (int i = 0; i < numTargetElems; ++i) { + delete targetElems[i]; + } + +} } diff --git a/src/INTERP_KERNEL/Interpolation3D.hxx b/src/INTERP_KERNEL/Interpolation3D.hxx index 069d9caaf..d21247438 100644 --- a/src/INTERP_KERNEL/Interpolation3D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D.hxx @@ -28,15 +28,22 @@ namespace MEDMEM public: // friend class INTERP_UTILS::MeshTestToolkit; - Interpolation3D(); - virtual ~Interpolation3D(); - - virtual IntersectionMatrix interpolateMeshes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh); - + Interpolation3D(){}; + virtual ~Interpolation3D(){}; + // Main function to interpolate triangular and quadratic meshes + template void interpolateMeshes(const MEDMEM::MESH& mesh1, + const MEDMEM::MESH& mesh2, + T& matrix); + //virtual IntersectionMatrix interpolateMeshes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh); + //MEDMEM::MEDMEM_Matrix* interpolate(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2); + void printInfo(int level) {}; + private: void calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, IntersectionMatrix& matrix); - + void calculateIntersectionVolumes(const MEDMEM::MESH& srcMesh, + const MEDMEM::MESH& targetMesh, + MEDMEM::MEDMEM_Matrix*& matrix); }; }; diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx new file mode 100644 index 000000000..dc4a9cef8 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D.txx @@ -0,0 +1,349 @@ +#include "Interpolation3D.hxx" +#include "MeshElement.hxx" +#include "TransformedTriangle.hxx" +#include "IntersectorTetra.hxx" +#include "IntersectorHexa.hxx" +#include "TargetIntersector.hxx" +#include "Log.hxx" + +using namespace INTERP_UTILS; +using namespace MEDMEM; +using namespace MED_EN; + +/// If defined, use recursion to traverse the binary search tree, else use the BBTree class +#define USE_RECURSIVE_BBOX_FILTER + +#ifdef USE_RECURSIVE_BBOX_FILTER +#include "MeshRegion.hxx" +#include "RegionNode.hxx" +#include + +#else // use BBTree class +#include "BBTree.H" + +#endif + +namespace MEDMEM { +/** + * \defgroup interpolation3D Interpolation3D + * \class Interpolation3D + * \brief Class used to calculate the volumes of intersection between the elements of two 3D meshes. + * + */ + +///** +// * Main interface method of the class, which verifies that the meshes are valid for the calculation, +// * and then returns the matrix of intersection volumes. +// * +// * @param srcMesh 3-dimensional source mesh +// * @param targetMesh 3-dimesional target mesh, containing only tetraedra +// * @return vector containing for each element i of the source mesh, a map giving for each element j +// * of the target mesh which i intersects, the volume of the intersection +// */ +//template Interpolation3D::interpolateMeshes( +// const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, T matrix&) { +// // we should do more sanity checking here - eliminating meshes that +// // cannot be treated +// +// calculateIntersectionVolumes(srcMesh, targetMesh, matrix); +// return matrix; +//} +// +//IntersectionMatrix Interpolation3D::interpolateMeshes( +// const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) { +// // it seems wasteful to make a copy here +// IntersectionMatrix matrix; +// +// // we should do more sanity checking here - eliminating meshes that +// // cannot be treated +// +// calculateIntersectionVolumes(srcMesh, targetMesh, matrix); +// return matrix; +//} + + +/** + * 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 by an object of type Intersector3D 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. + * 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 (#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 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param matrix vector of maps in which the result is stored + * + */ +template void Interpolation3D::interpolateMeshes( + const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh, T& matrix) { + + // create MeshElement objects corresponding to each element of the two meshes + const int numSrcElems = srcMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + const int numTargetElems = targetMesh.getNumberOfElements(MED_CELL, + MED_ALL_ELEMENTS); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + MeshElement* srcElems[numSrcElems]; + MeshElement* targetElems[numTargetElems]; + + std::map indices; + + for (int i = 0; i < numSrcElems; ++i) { + //const medGeometryElement type = srcMesh.getElementType(MED_CELL, i + 1); + srcElems[i] = new MeshElement(i + 1, srcMesh); + } + + for (int i = 0; i < numTargetElems; ++i) { + //const medGeometryElement type = targetMesh.getElementType(MED_CELL, i + 1); + targetElems[i] = new MeshElement(i + 1, targetMesh); + } + + // create empty maps for all source elements + matrix.resize(numTargetElems); + +#ifdef USE_RECURSIVE_BBOX_FILTER + + /* + * Performs a depth-first search over srcMesh, using bounding boxes to recursively eliminate the elements of targetMesh + * which cannot intersect smaller and smaller regions of srcMesh. At each level, each region is divided in two, forming + * a binary search tree with leaves consisting of only one element of the source mesh together with the elements of the + * target mesh that can intersect it. The recursion is implemented with a stack of RegionNodes, each one containing a + * source region and a target region. Each region has an associated bounding box and a vector of pointers to the elements + * that belong to it. Each MeshElement contains a bounding box and the global number of the corresponding element in the mesh. + */ + + // create initial RegionNode and fill up its source region with all the source mesh elements and + // its target region with all the target mesh elements whose bounding box + // intersects that of the source region + + RegionNode* firstNode = new RegionNode(); + + MeshRegion& srcRegion = firstNode->getSrcRegion(); + + for (int i = 0; i < numSrcElems; ++i) { + srcRegion.addElement(srcElems[i], srcMesh); + } + + MeshRegion& targetRegion = firstNode->getTargetRegion(); + + for (int i = 0; i < numTargetElems; ++i) { + if (!srcRegion.isDisjointWithElementBoundingBox( *(targetElems[i]))) { + targetRegion.addElement(targetElems[i], targetMesh); + } + } + + // Using a stack, descend recursively, creating at each step two new RegionNodes having as source region the left and + // right part of the source region of the current node (created using MeshRegion::split()) and as target region all the + // elements of the target mesh whose bounding box intersects the corresponding part + // Continue until the source region contains only one element, at which point the intersection volumes are + // calculated with all the remaining target mesh elements and stored in the matrix if they are non-zero. + + stack nodes; + nodes.push(firstNode); + + while (!nodes.empty()) { + RegionNode* currNode = nodes.top(); + nodes.pop();LOG(4, "Popping node "); + + if (currNode->getTargetRegion().getNumberOfElements() == 1) { + // calculate volumes + LOG(4, " - One element"); + + MeshElement* targetElement = *(currNode->getTargetRegion().getBeginElements()); + + // NB : srcElement indices are from 0 .. numSrcElements - 1 + // targetElement indicies from 1 .. numTargetElements + // maybe this is not ideal ... + const int targetIdx = targetElement->getIndex(); + + TargetIntersector* intersector; + + switch (targetElement->getType()) { + case MED_TETRA4: + intersector = new IntersectorTetra(srcMesh, targetMesh, targetIdx); + break; + + case MED_HEXA8: + intersector = new IntersectorHexa(srcMesh, targetMesh, targetIdx); + break; + + default: + assert(false); + } + + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { + + const int srcIdx = (*iter)->getIndex(); + const double vol = intersector->intersectSourceCell(srcIdx); + + if (vol != 0.0) { + matrix[targetIdx-1].insert(make_pair(srcIdx,vol)); + //matrix->setIJ(targetIdx,srcIdx,vol); + //matrix[targetIdx - 1][srcIdx] = vol; + //LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix->getIJ(srcIdx,targetIdx)); + + //LOG(2, "Result : V (" << srcIdx << ", " << targetIdx << ") = " << matrix[srcIdx - 1][targetIdx]); + + } + } + + delete intersector; + + } else // recursion + { + + LOG(4, " - Recursion"); + + RegionNode* leftNode = new RegionNode(); + RegionNode* rightNode = new RegionNode(); + + // split current source region + //} decide on axis + static BoundingBox::BoxCoord axis = BoundingBox::XMAX; + + currNode->getTargetRegion().split(leftNode->getTargetRegion(), + rightNode->getTargetRegion(), axis, targetMesh); + + LOG(5, "After split, left target region has " << leftNode->getTargetRegion().getNumberOfElements() + << " elements and right target region has " << rightNode->getTargetRegion().getNumberOfElements() + << " elements"); + + // ugly hack to avoid problem with enum which does not start at 0 + // I guess I ought to implement ++ for it instead ... + // Anyway, it basically chooses the next axis, cyclically + axis + = (axis != BoundingBox::ZMAX) ? static_cast(axis + + 1) + : BoundingBox::XMAX; + + // add source elements of current node that overlap the target regions of the new nodes + LOG(5, " -- Adding source elements"); + int numLeftElements = 0; + int numRightElements = 0; + for (vector::const_iterator iter = + currNode->getSrcRegion().getBeginElements() ; iter + != currNode->getSrcRegion().getEndElements() ; ++iter) { + LOG(6, " --- New target node"); + + if (!leftNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + leftNode->getSrcRegion().addElement(*iter, srcMesh); + ++numLeftElements; + } + + if (!rightNode->getTargetRegion().isDisjointWithElementBoundingBox(**iter)) { + rightNode->getSrcRegion().addElement(*iter, srcMesh); + ++numRightElements; + } + + } + + LOG(5, "Left src region has " << numLeftElements << " elements and right src region has " + << numRightElements << " elements"); + + // push new nodes on stack + if (numLeftElements != 0) { + nodes.push(leftNode); + } else { + delete leftNode; + } + + if (numRightElements != 0) { + nodes.push(rightNode); + } else { + delete rightNode; + } + } + + // all nodes are deleted here + delete currNode; + + LOG(4, "Next iteration. Nodes left : " << nodes.size()); + } + +#else // Use BBTree + // create BBTree structure + // - get bounding boxes + double bboxes[6 * numSrcElems]; + int srcElemIdx[numSrcElems]; + for(int i = 0; i < numSrcElems; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + + // source indices have to begin with zero for BBox, I think + srcElemIdx[i] = srcElems[i]->getIndex() - 1; + } + + BBTree<3> tree(bboxes, numSrcElems,srcElemIdx, 0); + + // for each target element, get source elements with which to calculate intersection + // - calculate intersection by calling intersectCells + for(int i = 0; i < numTargetElems; ++i) + { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + + vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + // create intersector + IntersectorTetra intersector(srcMesh, targetMesh, targetIdx); + + for(vector::const_iterator iter = intersectElems.begin(); iter != intersectElems.end(); ++iter) + { + + const int srcIdx = *iter + 1; + const double vol = intersector.intersectSourceCell(srcIdx); + + if(vol != 0.0) + { + matrix[targetIdx-1].insert(make_pair(srcIdx,vol)); + //matrix->setIJ(targetIdx,srcIdx,vol); + //matrix[targetIdx - 1].insert(make_pair(srcIdx, vol)); + } + + } + } + +#endif + + // free allocated memory + for (int i = 0; i < numSrcElems; ++i) { + delete srcElems[i]; + } + for (int i = 0; i < numTargetElems; ++i) { + delete targetElems[i]; + } + +} + + +} diff --git a/src/INTERP_KERNEL/Interpolation3DSurf.cxx b/src/INTERP_KERNEL/Interpolation3DSurf.cxx index 8994a8157..8b6f1c3b6 100644 --- a/src/INTERP_KERNEL/Interpolation3DSurf.cxx +++ b/src/INTERP_KERNEL/Interpolation3DSurf.cxx @@ -1,246 +1,57 @@ #include "MEDMEM_Mesh.hxx" #include "Interpolation3DSurf.hxx" -#include "PlanarIntersector.hxx" -#include "PlanarIntersector.H" -#include "TriangulationIntersector.hxx" -#include "TriangulationIntersector.cxx" -#include "ConvexIntersector.hxx" -#include "ConvexIntersector.cxx" -#include "GenericIntersector.hxx" -#include "GenericIntersector.cxx" #include "BBTree.H" #include - - using namespace MED_EN; -namespace MEDMEM -{ +namespace MEDMEM { /** -\defgroup interpolation3DSurf Interpolation3DSurf + \defgroup interpolation3DSurf Interpolation3DSurf \class Interpolation3DSurf -\brief Class used to compute the coefficients of the interpolation matrix between -two local surfacic meshes in three dimensions. Meshes can contain mixed triangular and quadrangular elements. -*/ - /** - * Constructor - * - */ - Interpolation3DSurf::Interpolation3DSurf() - { - _Intersection_type= MEDMEM::Triangulation; - _MedianPlane=0.5; - _Do_rotate=true; - _Precision=1.0E-12; - _DimCaracteristic=1; - _Surf3DAdjustmentEps=0.0001; - _PrintLevel=0; - } - /** - \brief Function used to set the options for the intersection calculation -\details The following options can be modified: + \brief Class used to compute the coefficients of the interpolation matrix between + two local surfacic meshes in three dimensions. Meshes can contain mixed triangular and quadrangular elements. + */ +/** + * Constructor + * + */ +Interpolation3DSurf::Interpolation3DSurf() { + _Intersection_type= MEDMEM::Triangulation; + _MedianPlane=0.5; + _Do_rotate=true; + _Precision=1.0E-12; + _DimCaracteristic=1; + _Surf3DAdjustmentEps=0.0001; + _PrintLevel=0; +} +/** + \brief Function used to set the options for the intersection calculation + \details The following options can be modified: -# Intersection_type: the type of algorithm to be used in the computation of the cell-cell intersections. - - Values: Triangle, Convex. - - Default: Triangle. + - Values: Triangle, Convex. + - Default: Triangle. -# MedianPlane: Position of the median plane where both cells will be projected - - Values: between 0 and 1. - - Default: 0.5. + - Values: between 0 and 1. + - Default: 0.5. -# DoRotate: rotate the coordinate system such that the target cell is in the Oxy plane. - - Values: true (necessarilly if Intersection_type=Triangle), false. - - Default: true (as default Intersection_type=Triangle) + - Values: true (necessarilly if Intersection_type=Triangle), false. + - Default: true (as default Intersection_type=Triangle) -# Precision: Level of precision of the computations is precision times the characteristic size of the mesh. - - Values: positive real number. - - Default: 1.0E-12. + - Values: positive real number. + - Default: 1.0E-12. -# PrintLevel: Level of verboseness during the computations. - - Values: interger between 0 and 3. - - Default: 0. + - Values: interger between 0 and 3. + - Default: 0. */ - void Interpolation3DSurf::setOptions(double Precision, int PrintLevel, double MedianPlane, IntersectionType intersection_type, bool do_rotate) - { - _Intersection_type=intersection_type; - _MedianPlane=MedianPlane; - _Do_rotate=do_rotate; - _Precision=Precision; - _PrintLevel=PrintLevel; - } - - - /** \brief Main function to interpolate triangular or quadrangular meshes. - \details The algorithm proceeds 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 by an object of type Intersector3Dsurf 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. - * 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 (#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 myMesh_S 3D surface source mesh - * @param myMesh_P 3D surface target mesh - * @return vector containing for each element i of the source mesh, a map giving for each element j - * of the target mesh which i intersects, the area of the intersection - * - */ - vector > Interpolation3DSurf::interpolateMeshes(const MEDMEM::MESH& myMesh_S, - const MEDMEM::MESH& myMesh_P) - { - long global_start =clock(); - /***********************************************************/ - /* Check both meshes are made of triangles and quadrangles */ - /***********************************************************/ - - int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL); - int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL); - if ( NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2) - cout<<"Interpolation3DSurf::interpolateMeshes: one of the two meshes contains more than two types of cells"< > BoxS=myMesh_S.getBoundingBox(); - vector > BoxP=myMesh_P.getBoundingBox(); - double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) - +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]) - +(BoxS[1][2]-BoxS[0][2])*(BoxS[1][2]-BoxS[0][2])); - double DimCaracteristic_S=diagonalS/nbMaille_S; - double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) - +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]) - +(BoxP[1][2]-BoxP[0][2])*(BoxP[1][2]-BoxP[0][2])); - double DimCaracteristic_P=diagonalP/nbMaille_P; - - _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); - if (_PrintLevel>=1) - { - cout<<" - Characteristic size of the source mesh : "<< DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, - _MedianPlane, _PrintLevel); - break; - case MEDMEM::Convex: - intersector=new ConvexIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, - _MedianPlane, _Do_rotate, _PrintLevel); - break; - case MEDMEM::Generic: - intersector=new GenericIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, - _MedianPlane, _Do_rotate, _PrintLevel); - break; - } - - /****************************************************************/ - /* Create a search tree based on the bounding boxes */ - /* Instanciate the intersector and initialise the result vector */ - /****************************************************************/ - - long start_filtering=clock(); - - vector bbox; - intersector->createBoundingBoxes<3>(myMesh_S,bbox); // create the bounding boxes - intersector->adjustBoundingBoxes<3>(bbox, _Surf3DAdjustmentEps); // expand them so that no intersection is missed - BBTree<3> my_tree(&bbox[0], 0, 0,nbMaille_S );//creating the search structure - - long end_filtering=clock(); - - map MAP_init; - vector > result_areas(nbMaille_P,MAP_init);//on initialise. - - /****************************************************/ - /* Loop on the target cells - core of the algorithm */ - /****************************************************/ - const MED_EN::medGeometryElement* types = myMesh_P.getTypes(MED_EN::MED_CELL); - int type_max_index=0;//maximum cell number for a given type - int type_min_index=0;//minimum cell number for a given type - int i_P=0;//global index of cell - - long start_intersection=clock(); - - for (int itype = 0; itype intersecting_elems; - double bb[6]; - intersector->getElemBB<3>(bb,myMesh_P,i_P+1,nb_nodesP); - my_tree.getIntersectingElems(bb, intersecting_elems); - int nb_intersecting_elems = intersecting_elems.size(); - //browsing all the i_S (from mesh S) elems that can - //intersect elem i_P (from mesh P) - for (int ielem=0; ielemintersectCells(i_P+1,i_S+1,nb_nodesP,nb_nodesS); - (result_areas[i_P]).insert(make_pair(i_S+1,surf)); - //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); - } - intersecting_elems.clear(); - } - type_min_index = type_max_index; - } - delete intersector; - - /***********************************/ - /* DEBUG prints */ - /***********************************/ - - if (_PrintLevel >=1) - { - long end_intersection=clock(); - if (_PrintLevel >=2) - { - cout<::iterator surface; - total_interm=0.0; - for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++) - { - cout<<" ("< + +using namespace MED_EN; + +namespace MEDMEM { +/** + \defgroup interpolation3DSurf Interpolation3DSurf + + \class Interpolation3DSurf + \brief Class used to compute the coefficients of the interpolation matrix between + two local surfacic meshes in three dimensions. Meshes can contain mixed triangular and quadrangular elements. + */ + + + +/** \brief Main function to interpolate triangular or quadrangular meshes. + \details The algorithm proceeds 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 by an object of type Intersector3Dsurf 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. + * 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 (#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 myMesh_S 3D surface source mesh + * @param myMesh_P 3D surface target mesh + * @return vector containing for each element i of the source mesh, a map giving for each element j + * of the target mesh which i intersects, the area of the intersection + * + */ +template void Interpolation3DSurf::interpolateMeshes( + const MEDMEM::MESH& myMesh_S, const MEDMEM::MESH& myMesh_P, T& matrix) { + long global_start =clock(); + /***********************************************************/ + /* Check both meshes are made of triangles and quadrangles */ + /***********************************************************/ + + int NumberOfCellsTypes_S = myMesh_S.getNumberOfTypes(MED_EN::MED_CELL); + int NumberOfCellsTypes_P = myMesh_P.getNumberOfTypes(MED_EN::MED_CELL); + if (NumberOfCellsTypes_S > 2 || NumberOfCellsTypes_P >2) + cout + <<"Interpolation3DSurf::interpolateMeshes: one of the two meshes contains more than two types of cells" + < > BoxS=myMesh_S.getBoundingBox(); + vector > BoxP=myMesh_P.getBoundingBox(); + double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) + +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]) +(BoxS[1][2] + -BoxS[0][2])*(BoxS[1][2]-BoxS[0][2])); + double DimCaracteristic_S=diagonalS/nbMaille_S; + double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) + +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]) +(BoxP[1][2] + -BoxP[0][2])*(BoxP[1][2]-BoxP[0][2])); + double DimCaracteristic_P=diagonalP/nbMaille_P; + + _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); + if (_PrintLevel>=1) { + cout<<" - Characteristic size of the source mesh : " + << DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, + _MedianPlane, _PrintLevel); + break; + case MEDMEM::Convex: + intersector=new ConvexIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, + _MedianPlane, _Do_rotate, _PrintLevel); + break; + case MEDMEM::Generic: + intersector=new GenericIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, + _MedianPlane, _Do_rotate, _PrintLevel); + break; + } + + /****************************************************************/ + /* Create a search tree based on the bounding boxes */ + /* Instanciate the intersector and initialise the result vector */ + /****************************************************************/ + + long start_filtering=clock(); + + vector bbox; + intersector->createBoundingBoxes<3>(myMesh_S, bbox); // create the bounding boxes + intersector->adjustBoundingBoxes<3>(bbox, _Surf3DAdjustmentEps); // expand them so that no intersection is missed + BBTree<3> my_tree(&bbox[0], 0, 0, nbMaille_S);//creating the search structure + + long end_filtering=clock(); + + map MAP_init; + //vector > result_areas(nbMaille_P, MAP_init);//on initialise. + matrix.resize(nbMaille_P); + /****************************************************/ + /* Loop on the target cells - core of the algorithm */ + /****************************************************/ + const MED_EN::medGeometryElement* types = + myMesh_P.getTypes(MED_EN::MED_CELL); + int type_max_index=0;//maximum cell number for a given type + int type_min_index=0;//minimum cell number for a given type + int i_P=0;//global index of cell + + long start_intersection=clock(); + + for (int itype = 0; itype intersecting_elems; + double bb[6]; + intersector->getElemBB<3>(bb, myMesh_P, i_P+1, nb_nodesP); + my_tree.getIntersectingElems(bb, intersecting_elems); + int nb_intersecting_elems = intersecting_elems.size(); + //browsing all the i_S (from mesh S) elems that can + //intersect elem i_P (from mesh P) + for (int ielem=0; ielemintersectCells(i_P+1, i_S+1, + nb_nodesP, nb_nodesS); + (matrix[i_P]).insert(make_pair(i_S+1, surf)); + //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); + } + intersecting_elems.clear(); + } + type_min_index = type_max_index; + } + delete intersector; + +// /***********************************/ +// /* DEBUG prints */ +// /***********************************/ +// +// if (_PrintLevel >=1) { +// long end_intersection=clock(); +// if (_PrintLevel >=2) { +// cout<::iterator surface; +// total_interm=0.0; +// for (surface=result_areas[i].begin(); surface +// !=result_areas[i].end(); surface++) { +// cout<<" ("< > BoxS=myMesh_S.getBoundingBox(); +// vector > BoxP=myMesh_P.getBoundingBox(); +// double diagonalS=sqrt((BoxS[1][0]-BoxS[0][0])*(BoxS[1][0]-BoxS[0][0]) +// +(BoxS[1][1]-BoxS[0][1])*(BoxS[1][1]-BoxS[0][1]) +(BoxS[1][2] +// -BoxS[0][2])*(BoxS[1][2]-BoxS[0][2])); +// double DimCaracteristic_S=diagonalS/nbMaille_S; +// double diagonalP=sqrt((BoxP[1][0]-BoxP[0][0])*(BoxP[1][0]-BoxP[0][0]) +// +(BoxP[1][1]-BoxP[0][1])*(BoxP[1][1]-BoxP[0][1]) +(BoxP[1][2] +// -BoxP[0][2])*(BoxP[1][2]-BoxP[0][2])); +// double DimCaracteristic_P=diagonalP/nbMaille_P; +// +// _DimCaracteristic=min(DimCaracteristic_S, DimCaracteristic_P); +// if (_PrintLevel>=1) { +// cout<<" - Characteristic size of the source mesh : " +// << DimCaracteristic_S<(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// _MedianPlane, _PrintLevel); +// break; +// case MEDMEM::Convex: +// intersector=new ConvexIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic, _Precision, +// _MedianPlane, _Do_rotate, _PrintLevel); +// break; +// case MEDMEM::Generic: +// intersector=new GenericIntersector<3>(myMesh_P,myMesh_S, _DimCaracteristic,_Precision, +// _MedianPlane, _Do_rotate, _PrintLevel); +// break; +// } +// +// /****************************************************************/ +// /* Create a search tree based on the bounding boxes */ +// /* Instanciate the intersector and initialise the result vector */ +// /****************************************************************/ +// +// long start_filtering=clock(); +// +// vector bbox; +// intersector->createBoundingBoxes<3>(myMesh_S, bbox); // create the bounding boxes +// intersector->adjustBoundingBoxes<3>(bbox, _Surf3DAdjustmentEps); // expand them so that no intersection is missed +// BBTree<3> my_tree(&bbox[0], 0, 0, nbMaille_S);//creating the search structure +// +// long end_filtering=clock(); +// +// map MAP_init; +// +// MEDMEM_Matrix* matrix = new MEDMEM_Matrix (nbMaille_P); +// //vector > result_areas(nbMaille_P, MAP_init);//on initialise. +// +// /****************************************************/ +// /* Loop on the target cells - core of the algorithm */ +// /****************************************************/ +// const MED_EN::medGeometryElement* types = +// myMesh_P.getTypes(MED_EN::MED_CELL); +// int type_max_index=0;//maximum cell number for a given type +// int type_min_index=0;//minimum cell number for a given type +// int i_P=0;//global index of cell +// +// long start_intersection=clock(); +// +// for (int itype = 0; itype intersecting_elems; +// double bb[6]; +// intersector->getElemBB<3>(bb, myMesh_P, i_P+1, nb_nodesP); +// my_tree.getIntersectingElems(bb, intersecting_elems); +// int nb_intersecting_elems = intersecting_elems.size(); +// //browsing all the i_S (from mesh S) elems that can +// //intersect elem i_P (from mesh P) +// for (int ielem=0; ielemintersectCells(i_P+1, i_S+1, +// nb_nodesP, nb_nodesS); +// //(result_areas[i_P]).insert(make_pair(i_S+1, surf)); +// matrix->setIJ(i_P+1,i_S+1,surf); +// //rempli_vect_d_intersect(i_P+1,i_S+1,MaP,MaS,result_areas); +// } +// intersecting_elems.clear(); +// } +// type_min_index = type_max_index; +// } +// delete intersector; +// +// return matrix; +//} +} +; diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index ebd5310c5..29041ce75 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -90,11 +90,11 @@ namespace INTERP_UTILS /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ - /* fonction qui calcul le déterminant */ + /* fonction qui calcul le d�terminant */ /* de deux vecteur(cf doc CGAL). */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ - //fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2 + //fonction qui calcul le d�terminant des vecteurs: P3P1 et P3P2 //(cf doc CGAL). inline double mon_determinant(const double* P_1, @@ -146,8 +146,8 @@ namespace INTERP_UTILS /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ - /* calcul les coordonnées du barycentre d'un polygone */ - /* le vecteur en entrée est constitué des coordonnées */ + /* calcul les coordonn�es du barycentre d'un polygone */ + /* le vecteur en entr�e est constitu� des coordonn�es */ /* des sommets du polygone */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -180,7 +180,7 @@ namespace INTERP_UTILS { double Surface=0; - for (int i=0; i<(Poly.size())/2-2; i++) + for (uint i=0; i<(Poly.size())/2-2; i++) { double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] ); Surface=Surface + Surf ; @@ -215,7 +215,7 @@ namespace INTERP_UTILS /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ - /*fonction pour vérifier qu'un point n'a pas déja été considérer dans */ + /*fonction pour v�rifier qu'un point n'a pas d�ja �t� consid�rer dans */ /* le vecteur et le rajouter au vecteur sinon. */ /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ @@ -238,9 +238,9 @@ namespace INTERP_UTILS // V.push_back(P[0]); // V.push_back(P[1]); - int taille=V.size(); + uint taille=V.size(); bool isPresent=false; - for(int i=0;i tab; - int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n°des noeuds commencent à 1. + int I_1=MaS._ReversNConnectIndex[N_1-1]-1;// Dans MED le n�des noeuds commencent � 1. int I_2=MaS._ReversNConnectIndex[N_2-1]-1; int I_3=MaS._ReversNConnectIndex[N_3-1]-1; int F_1=MaS._ReversNConnectIndex[N_1]-1; @@ -402,7 +402,7 @@ namespace INTERP_UTILS MaS._ReversNConnect+F_3,ib_3); - for(int i=0;i<(V_12).size();i++) + for(uint i=0;i<(V_12).size();i++) { if(V_12[i]!=i_S) { @@ -411,7 +411,7 @@ namespace INTERP_UTILS } } - for(int i=0;i& FractSurf,MEDMEM::MESH* myMesh_P, MEDMEM::MESH* myMesh_S,string NameMesh_P,string NameMesh_S) { - //On sauvegarde le champ crée sur la maillage P. + //On sauvegarde le champ cr�e sur la maillage P. string FileName="FracSurf"+NameMesh_S+"_to_"+NameMesh_P; int id1=(*myMesh_P).addDriver(MED_DRIVER,FileName,"MP"); int id2=(*myMesh_S).addDriver(MED_DRIVER,FileName,"MS"); diff --git a/src/INTERP_KERNEL/IntersectorTetra.cxx b/src/INTERP_KERNEL/IntersectorTetra.cxx index 2df69a7d5..6bcb71920 100644 --- a/src/INTERP_KERNEL/IntersectorTetra.cxx +++ b/src/INTERP_KERNEL/IntersectorTetra.cxx @@ -202,10 +202,10 @@ namespace INTERP_UTILS // calculate the triangles if needed // local nodes 1, 2, 3 - TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[2]); + TriangleFaceKey key1 = TriangleFaceKey(faceNodes[0], faceNodes[1], faceNodes[3]); if(_volumes.find(key1) == _volumes.end()) { - TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[2]]); + TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[1]], _nodes[faceNodes[3]]); calculateVolume(tri, key1); totalVolume += _volumes[key1]; } else { @@ -214,10 +214,10 @@ namespace INTERP_UTILS } // local nodes 1, 3, 4 - TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); + TriangleFaceKey key2 = TriangleFaceKey(faceNodes[1], faceNodes[2], faceNodes[3]); if(_volumes.find(key2) == _volumes.end()) { - TransformedTriangle tri(_nodes[faceNodes[0]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]); + TransformedTriangle tri(_nodes[faceNodes[1]], _nodes[faceNodes[2]], _nodes[faceNodes[3]]); calculateVolume(tri, key2); totalVolume += _volumes[key2]; } @@ -225,8 +225,9 @@ namespace INTERP_UTILS { // count negative as face has reversed orientation totalVolume -= _volumes[key2]; - } + } } + break; default: diff --git a/src/INTERP_KERNEL/Makefile.in b/src/INTERP_KERNEL/Makefile.in index fd38c3a8a..3ecceaf35 100644 --- a/src/INTERP_KERNEL/Makefile.in +++ b/src/INTERP_KERNEL/Makefile.in @@ -66,7 +66,8 @@ TriangulationIntersector.cxx\ ConvexIntersector.hxx\ ConvexIntersector.cxx\ GenericIntersector.hxx\ -GenericIntersector.cxx +GenericIntersector.cxx\ +Remapper.hxx # Libraries targets @@ -82,12 +83,12 @@ MeshRegion.cxx\ RegionNode.cxx\ BoundingBox.cxx\ TetraAffineTransform.cxx\ -Interpolation3D.cxx\ Interpolation3DSurf.cxx\ Interpolation2D.cxx\ IntersectorTetra.cxx\ IntersectorHexa.cxx\ -PlanarIntersector.cxx +PlanarIntersector.cxx\ +Remapper.cxx # Executables targets BIN = diff --git a/src/INTERP_KERNEL/Remapper.cxx b/src/INTERP_KERNEL/Remapper.cxx new file mode 100644 index 000000000..b926f968b --- /dev/null +++ b/src/INTERP_KERNEL/Remapper.cxx @@ -0,0 +1,113 @@ +#include "Remapper.hxx" +#include "MEDMEM_Exception.hxx" +#include "Interpolation.hxx" +#include "Interpolation2D.hxx" +#include "Interpolation3D.hxx" +#include "Interpolation3DSurf.hxx" +#include "Interpolation2D.txx" +#include "Interpolation3D.txx" +#include "Interpolation3DSurf.txx" + +namespace INTERP_KERNEL +{ + +Remapper::Remapper():_matrix(0) +{ +} + +Remapper::~Remapper() +{ + delete _matrix; +} +/*! This method computes the intersection matrix between + * source \a mesh_source and \a mesh_target. It is a preliminary step + * that is necessary before calling the \c transfer() method. + * The method analyses the dimensions of the meshes and checks for compatibility. + * + */ +void Remapper::prepare(const MEDMEM::MESH& mesh_source, const MEDMEM::MESH& mesh_target) +{ + int sm_spacedim = mesh_source.getSpaceDimension(); + int tm_spacedim = mesh_target.getSpaceDimension(); + int sm_meshdim = mesh_source.getMeshDimension(); + int tm_meshdim = mesh_target.getMeshDimension(); + + if (tm_spacedim!=sm_spacedim || tm_meshdim!=sm_meshdim) + throw MEDEXCEPTION("incompatible mesh and/or space dimensions in meshes"); + + _matrix= new MEDMEM::MEDMEM_Matrix; + if ((sm_spacedim==2)&&(sm_meshdim==2)) + { + Interpolation2D interpolation; + interpolation.interpolateMeshes(mesh_source,mesh_target,*_matrix); + } + else if ((sm_spacedim==3)&&(sm_meshdim==3)) + { + Interpolation3D interpolation; + interpolation.interpolateMeshes(mesh_source,mesh_target,*_matrix); + } + else if ((sm_spacedim==3)&&(sm_meshdim==2)) + { + Interpolation3DSurf interpolation; + interpolation.interpolateMeshes(mesh_source,mesh_target,*_matrix); + } + else + throw MEDEXCEPTION("no Interpolation exists for the given mesh and space dimensions"); + +} + +void Remapper::transfer(const MEDMEM::FIELD& field_source,MEDMEM::FIELD& field_target) +{ + int source_nbcomp=field_source.getNumberOfComponents(); + int target_nbcomp=field_target.getNumberOfComponents(); + if (source_nbcomp != target_nbcomp) + throw MEDMEM::MEDEXCEPTION("incoherent number of components for source and target fields"); + if (source_nbcomp>1) + throw MEDMEM::MEDEXCEPTION("interpolations with more than one component are not yet handled"); + MEDMEM::FIELD* target_volumes = getSupportVolumes(*field_target.getSupport()); + int nbelem_target=field_target.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS); + + double* value_target = const_cast (field_target.getValue()); + const double* value_source = field_source.getValue(); + + _matrix->multiply(value_source, value_target); + for (int i=0; i< nbelem_target; i++) + value_target[i]/=target_volumes->getValueIJ(i+1,1); + + delete target_volumes; +} + +void Remapper::setOptionDouble(const std::string& key, double value) +{ +} + +void Remapper::setOptionInt(const std::string& key, int value) +{ +} + +/*! +\brief returns the volumes of the cells underlying the field \a field + +For 2D geometries, the returned field contains the areas. +For 3D geometries, the returned field contains the volumes. + +\param field field on which cells the volumes are required +\return field containing the volumes +*/ + MEDMEM::FIELD* Remapper::getSupportVolumes(const MEDMEM::SUPPORT& support) + { + const MEDMEM::MESH* mesh=support.getMesh(); + int dim = mesh->getMeshDimension(); + switch (dim) + { + case 2: + return mesh->getArea(&support); + case 3: + return mesh->getVolume(&support); + default: + throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension"); + } +} + + +} diff --git a/src/INTERP_KERNEL/Remapper.hxx b/src/INTERP_KERNEL/Remapper.hxx new file mode 100644 index 000000000..026eab67a --- /dev/null +++ b/src/INTERP_KERNEL/Remapper.hxx @@ -0,0 +1,29 @@ +#ifndef REMAPPER_HXX_ +#define REMAPPER_HXX_ + +#include "MEDMEM_Matrix.hxx" +#include "MEDMEM_Mesh.hxx" +#include "MEDMEM_Support.hxx" +#include "MEDMEM_Field.hxx" + +namespace INTERP_KERNEL +{ + +class Remapper +{ +public: + Remapper(); + virtual ~Remapper(); + void prepare(const MEDMEM::MESH& mesh_source, const MEDMEM::MESH& mesh_target); + void transfer(const MEDMEM::FIELD& field_source, MEDMEM::FIELD& field_target); + void setOptionDouble(const std::string& key, double value); + void setOptionInt(const std::string& key, int value); +private : + MEDMEM::MEDMEM_Matrix* _matrix; + MEDMEM::FIELD* getSupportVolumes(const MEDMEM::SUPPORT& support); + +}; + +} + +#endif /*REMAPPER_HXX_*/ diff --git a/src/INTERP_KERNEL/Test/Makefile.in b/src/INTERP_KERNEL/Test/Makefile.in index e3b27d019..62c8d2c01 100644 --- a/src/INTERP_KERNEL/Test/Makefile.in +++ b/src/INTERP_KERNEL/Test/Makefile.in @@ -40,8 +40,9 @@ EXPORT_HEADERS = CppUnitTest.hxx \ Interpolation3DTestSuite.hxx \ SingleElementTetraTests.hxx \ MultiElementTetraTests.hxx \ - HexaTests.hxx\ - MeshTestToolkit.hxx + HexaTests.hxx \ + MeshTestToolkit.hxx \ + BBTreeTest.hxx @@ -59,7 +60,12 @@ TEST_PROGS = TestInterpKernel BIN = PerfTest -BIN_SRC = MeshTestToolkit.cxx CppUnitTest.cxx TransformedTriangleTest.cxx TransformedTriangleIntersectTest.cxx MeshTestToolkit.cxx +BIN_SRC = MeshTestToolkit.cxx \ + CppUnitTest.cxx \ + TransformedTriangleTest.cxx \ + TransformedTriangleIntersectTest.cxx \ + MeshTestToolkit.cxx \ + BBTreeTest.cxx BIN_CLIENT_IDL = @@ -76,16 +82,18 @@ CXXFLAGS+= $(CPPUNIT_INCLUDES) CPPFLAGS+= $(CPPUNIT_INCLUDES) # for log -CXXFLAGS+= -DLOG_LEVEL=3 -DOPTIMIZE -O2 -Wall -CPPFLAGS+= -DLOG_LEVEL=3 -DOPTIMIZE -O2 -Wall +CXXFLAGS+= -DLOG_LEVEL=3 -DOPTIMIZE -Wall +CPPFLAGS+= -DLOG_LEVEL=3 -DOPTIMIZE -Wall +#CXXFLAGS += -O2 +#CPPFLAGS += -O2 # for gcov #CXXFLAGS+=-fprofile-arcs -ftest-coverage #CPPFLAGS+=-fprofile-arcs -ftest-coverage #for gprof -CXXFLAGS+=-pg -CPPFLAGS+=-pg +#CXXFLAGS+=-pg +#CPPFLAGS+=-pg #LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. @@ -93,8 +101,8 @@ LDFLAGS+=$(MED2_LIBS) $(HDF5_LIBS) -lmed_V2_1 $(STDLIB) -lmedmem -lutil #LDFLAGS+= $(HDF5_LIBS) $(STDLIB) -lmedmem -lutil -lmed # for gcov -LDFLAGS+= -lgcov -LDFLAGS+= -pg +#LDFLAGS+= -lgcov +#LDFLAGS+= -pg #LDFLAGSFORBIN+=$(MED2_LIBS) $(HDF5_LIBS) # change motivated by the bug KERNEL4778. diff --git a/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx b/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx index 38eb3f11a..1b13e0479 100644 --- a/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx +++ b/src/INTERP_KERNEL/Test/MeshTestToolkit.cxx @@ -1,5 +1,6 @@ #include "MeshTestToolkit.hxx" #include "MEDMEM_Mesh.hxx" +#include "Interpolation3D.txx" #include #include @@ -328,7 +329,7 @@ namespace INTERP_TEST LOG(5, "Loading " << mesh2 << " from " << mesh2path); MESH tMesh(MED_DRIVER, dataDir+mesh2path, mesh2); - m = _interpolator->interpolateMeshes(sMesh, tMesh); + _interpolator->interpolateMeshes(sMesh, tMesh,m); // if reflexive, check volumes if(strcmp(mesh1path,mesh2path) == 0) diff --git a/src/INTERP_KERNEL/Test/PerfTest.cxx b/src/INTERP_KERNEL/Test/PerfTest.cxx index fd06dec37..391278b08 100644 --- a/src/INTERP_KERNEL/Test/PerfTest.cxx +++ b/src/INTERP_KERNEL/Test/PerfTest.cxx @@ -1,4 +1,5 @@ #include "Interpolation3D.hxx" +#include "Interpolation3D.txx" #include "MeshTestToolkit.hxx" #include "Log.hxx" #include "VectorUtils.hxx" @@ -63,7 +64,7 @@ namespace INTERP_TEST LOG(1, "Target mesh has " << numTargetElems << " elements"); - m = _interpolator->interpolateMeshes(sMesh, tMesh); + _interpolator->interpolateMeshes(sMesh, tMesh,m); std::pair eff = countNumberOfMatrixEntries(m); LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = " diff --git a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx index 6252407c0..a20a274bd 100644 --- a/src/INTERP_KERNEL/Test/TestInterpKernel.cxx +++ b/src/INTERP_KERNEL/Test/TestInterpKernel.cxx @@ -24,6 +24,7 @@ #include "MultiElementTetraTests.hxx" #include "SingleElementTetraTests.hxx" #include "HexaTests.hxx" +#include "BBTreeTest.hxx" using namespace INTERP_TEST; @@ -33,7 +34,7 @@ CPPUNIT_TEST_SUITE_REGISTRATION( MultiElementTetraTests ); CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementTetraTests ); CPPUNIT_TEST_SUITE_REGISTRATION( INTERP_TEST::TransformedTriangleIntersectTest ); CPPUNIT_TEST_SUITE_REGISTRATION( INTERP_TEST::TransformedTriangleTest ); - +CPPUNIT_TEST_SUITE_REGISTRATION( BBTreeTest); // --- generic Main program from KERNEL_SRC/src/Basics/Test #include "BasicMainTest.hxx" diff --git a/src/INTERP_KERNEL/TransformedTriangle_math.cxx b/src/INTERP_KERNEL/TransformedTriangle_math.cxx index c9928202e..f8f3035d5 100644 --- a/src/INTERP_KERNEL/TransformedTriangle_math.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle_math.cxx @@ -49,7 +49,7 @@ namespace INTERP_UTILS const long double TransformedTriangle::MULT_PREC_F = 4.0 * TransformedTriangle::MACH_EPS; /// Threshold for resetting double and triple products to zero; ( F / f in Grandy ) - const long double TransformedTriangle::THRESHOLD_F = 500.0; + const long double TransformedTriangle::THRESHOLD_F = 2000.0; /// Threshold for what is considered a small enough angle to warrant correction of triple products by Grandy, [57] const double TransformedTriangle::TRIPLE_PRODUCT_ANGLE_THRESHOLD = 0.1;