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)
}
_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);
}
bool intersects = true;
for (int idim=0; idim<dim; idim++)
{
- if (bb_ptr[idim*2]-bb[idim*2+1]>1e-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)
#include <map>
#include "MEDMEM_Mesh.hxx"
-
+#include "MEDMEM_Matrix.hxx"
/**
* \mainpage
* Status : documentation of 3D - part of intersection matrix calculation more or less complete
Interpolation(){}
virtual ~Interpolation(){}
- //interpolation of two triangular meshes.
- virtual std::vector<std::map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)=0;
-
+ virtual void printInfo(int level)=0;
};
};
#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<time.h>
using namespace MED_EN;
_PrintLevel=PrintLevel;
_Intersection_type=intersection_type;
}
+// std::vector<map<int, double> > Interpolation2D::interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)
+// {
+// std::vector<map<int,double> > matrix;
+// interpolateMeshes(mesh1,mesh2,matrix);
+// return matrix;
+// }
+//
+// MEDMEM::MEDMEM_Matrix<double,1> Interpolation2D::interpolate(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2)
+// {
+// MEDMEM::MEDMEM_Matrix<double,1> 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<map<int,double> > 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"<<endl;
- string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
- string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
-
- if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
- cout<<"Interpolation2D::intersectMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
- //VB
- delete[] Type_S;
- delete[] Type_P;
- int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
- int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
-
- /**************************************************/
- /* Search the characteristic size of the meshes */
- /**************************************************/
-
- vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
- vector<vector<double> > 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<<endl;
- cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
- cout << "Interpolation2D::computation of the intersections"<<endl;
- }
-
- PlanarIntersector* intersector;
-
- switch (_Intersection_type)
- {
- case MEDMEM::Triangulation:
- intersector=new TriangulationIntersector<2>(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<double> 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 <class T>
+// 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"<<endl;
+// string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+// string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+//
+// if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+// cout<<"Interpolation2D::intersectMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
+// //VB
+// delete[] Type_S;
+// delete[] Type_P;
+// int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+// int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+//
+// /**************************************************/
+// /* Search the characteristic size of the meshes */
+// /**************************************************/
+//
+// vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+// vector<vector<double> > 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<<endl;
+// cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
+// cout << "Interpolation2D::computation of the intersections"<<endl;
+// }
+//
+// PlanarIntersector* intersector;
+//
+// switch (_Intersection_type)
+// {
+// case MEDMEM::Triangulation:
+// intersector=new TriangulationIntersector<2>(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<double> 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<int,double> MAP_init;
+// vector<map<int,double> > 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<NumberOfCellsTypes_P; itype++)
+// {
+// int nb_nodesP=types[itype]%100;
+// int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+// type_max_index += nbelem_type;
+// for( i_P=type_min_index; i_P<type_max_index ; i_P++)
+// {
+// vector<int> 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; ielem<nb_intersecting_elems;ielem++)
+// {
+// //BBTree structure returns numbers between 0 and n-1
+// int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+// int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
+// double surf=intersector->intersectCells(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<<endl<<"Printing intersection areas:"<<endl<<endl;
+// cout<<"(source cell, target cell): intersection areas"<<endl;
+// double total=0.0;
+// double total_interm=0.0;
+// int nb_result_areas = result_areas.size();
+// for(int i=0; i< nb_result_areas;i++)
+// {
+// map<int,double>::iterator surface;
+// total_interm=0.0;
+// for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
+// {
+// cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+// total_interm +=(*surface).second;
+// }
+// cout<< " elem " << i+1 << " area= " << total_interm << endl;
+// total+=total_interm;
+// }
+// cout << "total area "<<total<<endl;
+// }
+// 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 result_areas;
+// }
+//
+// /** \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 sparse : it is a MEDMEM::Matrix<double,1> structure, i.e. a sparse matrix with 1-indexing
+// *
+//
+// * @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
+// *
+// */
+//
+// MEDMEM::MEDMEM_Matrix<double,1>* Interpolation2D::interpolate(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"<<endl;
+// string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+// string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+//
+// if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+// cout<<"Interpolation2D::intersectMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
+// //VB
+// delete[] Type_S;
+// delete[] Type_P;
+// int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+// int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+//
+// /**************************************************/
+// /* Search the characteristic size of the meshes */
+// /**************************************************/
+//
+// vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+// vector<vector<double> > 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<<endl;
+// cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
+// cout << "Interpolation2D::computation of the intersections"<<endl;
+// }
+//
+// PlanarIntersector* intersector;
+//
+// switch (_Intersection_type)
+// {
+// case MEDMEM::Triangulation:
+// intersector=new TriangulationIntersector<2>(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<double> 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<int,double> MAP_init;
+// MEDMEM_Matrix<double,1>* matrix=new MEDMEM_Matrix<double,1>(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<NumberOfCellsTypes_P; itype++)
+// {
+// int nb_nodesP=types[itype]%100;
+// int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+// type_max_index += nbelem_type;
+// for( i_P=type_min_index; i_P<type_max_index ; i_P++)
+// {
+// vector<int> 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; ielem<nb_intersecting_elems;ielem++)
+// {
+// //BBTree structure returns numbers between 0 and n-1
+// int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+// int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
+// double surf=intersector->intersectCells(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<int,double> MAP_init;
- vector<map<int,double> > 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<NumberOfCellsTypes_P; itype++)
- {
- int nb_nodesP=types[itype]%100;
- int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
- type_max_index += nbelem_type;
- for( i_P=type_min_index; i_P<type_max_index ; i_P++)
- {
- vector<int> 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; ielem<nb_intersecting_elems;ielem++)
- {
- //BBTree structure returns numbers between 0 and n-1
- int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
- int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
- double surf=intersector->intersectCells(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<<endl<<"Printing intersection areas:"<<endl<<endl;
- cout<<"(source cell, target cell): intersection areas"<<endl;
- double total=0.0;
- double total_interm=0.0;
- int nb_result_areas = result_areas.size();
- for(int i=0; i< nb_result_areas;i++)
- {
- map<int,double>::iterator surface;
- total_interm=0.0;
- for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
- {
- cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
- total_interm +=(*surface).second;
- }
- cout<< " elem " << i+1 << " area= " << total_interm << endl;
- total+=total_interm;
- }
- cout << "total area "<<total<<endl;
- }
- 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 result_areas;
- }
};
void setOptions(double Precision, int PrintLevel,
IntersectionType intersectionType);
- // Main function to interpolate triangular and qudadratic meshes
- std::vector<map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+ // Main function to interpolate triangular and quadratic meshes
+ template<class T> void interpolateMeshes(const MEDMEM::MESH& mesh1,
+ const MEDMEM::MESH& mesh2,
+ T& matrix);
+ void printInfo(int level) {};
private:
};
--- /dev/null
+#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"
+
+using namespace MED_EN;
+
+namespace MEDMEM {
+ /** \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 <class T>
+ 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"<<endl;
+ string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+ string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+
+ if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+ cout<<"Interpolation2D::intersectMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
+ //VB
+ delete[] Type_S;
+ delete[] Type_P;
+ int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+ int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
+
+ /**************************************************/
+ /* Search the characteristic size of the meshes */
+ /**************************************************/
+
+ vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+ vector<vector<double> > 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<<endl;
+ cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
+ cout << "Interpolation2D::computation of the intersections"<<endl;
+ }
+
+ PlanarIntersector* intersector;
+
+ switch (_Intersection_type)
+ {
+ case MEDMEM::Triangulation:
+ intersector=new TriangulationIntersector<2>(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<double> 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<NumberOfCellsTypes_P; itype++)
+ {
+ int nb_nodesP=types[itype]%100;
+ int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
+ type_max_index += nbelem_type;
+ for( i_P=type_min_index; i_P<type_max_index ; i_P++)
+ {
+ vector<int> 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; ielem<nb_intersecting_elems;ielem++)
+ {
+ //BBTree structure returns numbers between 0 and n-1
+ int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+ int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
+ double surf=intersector->intersectCells(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<<endl<<"Printing intersection areas:"<<endl<<endl;
+// cout<<"(source cell, target cell): intersection areas"<<endl;
+// double total=0.0;
+// double total_interm=0.0;
+// int nb_result_areas = result_areas.size();
+// for(int i=0; i< nb_result_areas;i++)
+// {
+// map<int,double>::iterator surface;
+// total_interm=0.0;
+// for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
+// {
+// cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+// total_interm +=(*surface).second;
+// }
+// cout<< " elem " << i+1 << " area= " << total_interm << endl;
+// total+=total_interm;
+// }
+// cout << "total area "<<total<<endl;
+// }
+ 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;
+ }
+
+ }
+
+};
#include "TargetIntersector.hxx"
#include "Log.hxx"
-
using namespace INTERP_UTILS;
using namespace MEDMEM;
using namespace MED_EN;
#include <stack>
#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<double,1>* Interpolation3D::interpolate(
+ const MEDMEM::MESH& srcMesh, const MEDMEM::MESH& targetMesh) {
+ // it seems wasteful to make a copy here
+ //IntersectionMatrix matrix
+ MEDMEM::MEDMEM_Matrix<double,1>* 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<MeshElement*, int> 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<double,1>*& 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<MeshElement*, int> 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<double,1>(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<RegionNode*> 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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<RegionNode*> 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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<int> intersectElems;
-
- tree.getIntersectingElems(targetBox, intersectElems);
-
- // create intersector
- IntersectorTetra intersector(srcMesh, targetMesh, targetIdx);
-
- for(vector<int>::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<int> intersectElems;
+
+ tree.getIntersectingElems(targetBox, intersectElems);
+
+ // create intersector
+ IntersectorTetra intersector(srcMesh, targetMesh, targetIdx);
+
+ for(vector<int>::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<MeshElement*, int> 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<RegionNode*> 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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<int> intersectElems;
+
+ tree.getIntersectingElems(targetBox, intersectElems);
+
+ // create intersector
+ IntersectorTetra intersector(srcMesh, targetMesh, targetIdx);
+
+ for(vector<int>::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];
+ }
+
+}
}
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<class T> 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<double,1>* 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<double,1>*& matrix);
};
};
--- /dev/null
+#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 <stack>
+
+#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 <class T> 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 <class T> 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<MeshElement*, int> 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<RegionNode*> 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<MeshElement*>::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<BoundingBox::BoxCoord>(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<MeshElement*>::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<int> intersectElems;
+
+ tree.getIntersectingElems(targetBox, intersectElems);
+
+ // create intersector
+ IntersectorTetra intersector(srcMesh, targetMesh, targetIdx);
+
+ for(vector<int>::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];
+ }
+
+}
+
+
+}
#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<time.h>
-
-
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<map<int,double> > 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"<<endl;
- string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
- string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
-
- if((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0] != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
- cout<<"Interpolation3DSurf::interpolateMeshes: one of the two meshes is neither linear triangular nor linear quadratic"<<endl;
- //VB
- delete[] Type_S;
- delete[] Type_P;
- int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
- int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,MED_EN::MED_ALL_ELEMENTS);
-
- /**************************************************/
- /* Search the characteristic size of the meshes */
- /**************************************************/
-
- vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
- vector<vector<double> > 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<<endl;
- cout<<" - Characteristic size of the target mesh: "<< DimCaracteristic_P<<endl;
- cout << "Interpolation3DSurf::computation of the intersections"<<endl;
- }
-
- PlanarIntersector* intersector;
-
- switch (_Intersection_type)
- {
- case MEDMEM::Triangulation:
- intersector=new TriangulationIntersector<3>(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<double> 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<int,double> MAP_init;
- vector<map<int,double> > 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<NumberOfCellsTypes_P; itype++)
- {
- int nb_nodesP=types[itype]%100;
- int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL, types[itype]);
- type_max_index += nbelem_type;
- for( i_P=type_min_index; i_P<type_max_index ; i_P++)
- {
- vector<int> 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; ielem<nb_intersecting_elems;ielem++)
- {
- //BBTree structure returns numbers between 0 and n-1
- int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
- int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL,i_S+1)%100;
- double surf=intersector->intersectCells(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<<endl<<"Printing intersection areas:"<<endl<<endl;
- cout<<"(source cell, target cell): intersection areas"<<endl;
- double total=0.0;
- double total_interm=0.0;
- int nb_result_areas = result_areas.size();
- for(int i=0; i< nb_result_areas;i++)
- {
- map<int,double>::iterator surface;
- total_interm=0.0;
- for( surface=result_areas[i].begin();surface!=result_areas[i].end();surface++)
- {
- cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
- total_interm +=(*surface).second;
- }
- cout<< " elem " << i+1 << " area= " << total_interm << endl;
- total+=total_interm;
- }
- cout << "total area "<<total<<endl;
- }
- cout<< "Filtering time= " << end_filtering-start_filtering<< endl;
- cout<< "Intersection time= " << end_intersection-start_intersection<< endl;
- long global_end =clock();
- cout<< "Global time= " << global_end - global_start << endl;
- }
-
- return result_areas;
- }
-};
+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;
+}
+}
+;
void setOptions(double Precision, int PrintLevel, double MedianPlane,
IntersectionType intersectionType, bool do_rotate);
- // Main function to interpolate triangular and qudadratic meshes
- std::vector<map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
-
+ // Main function to interpolate triangular and quadratic meshes
+ template<class T> void interpolateMeshes(const MEDMEM::MESH& mesh1,
+ const MEDMEM::MESH& mesh2,
+ T& matrix);
+// // Main function to interpolate triangular and qudadratic meshes
+// std::vector<map<int, double> > interpolateMeshes(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+// MEDMEM::MEDMEM_Matrix<double,1>* interpolate(const MEDMEM::MESH& mesh1, const MEDMEM::MESH& mesh2);
+ void printInfo(int level) {};
+
private:
};
--- /dev/null
+#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<time.h>
+
+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 <class T> 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"
+ <<endl;
+ string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+ string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+
+ if ((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0]
+ != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+ cout
+ <<"Interpolation3DSurf::interpolateMeshes: one of the two meshes is neither linear triangular nor linear quadratic"
+ <<endl;
+ //VB
+ delete[] Type_S;
+ delete[] Type_P;
+ int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,
+ MED_EN::MED_ALL_ELEMENTS);
+ int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,
+ MED_EN::MED_ALL_ELEMENTS);
+
+ /**************************************************/
+ /* Search the characteristic size of the meshes */
+ /**************************************************/
+
+ vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+ vector<vector<double> > 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<<endl;
+ cout<<" - Characteristic size of the target mesh: "
+ << DimCaracteristic_P<<endl;
+ cout << "Interpolation3DSurf::computation of the intersections"<<endl;
+ }
+
+ PlanarIntersector* intersector;
+
+ switch (_Intersection_type) {
+ case MEDMEM::Triangulation:
+ intersector=new TriangulationIntersector<3>(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<double> 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<int,double> MAP_init;
+ //vector<map<int,double> > 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<NumberOfCellsTypes_P; itype++) {
+ int nb_nodesP=types[itype]%100;
+ int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL,
+ types[itype]);
+ type_max_index += nbelem_type;
+ for (i_P=type_min_index; i_P<type_max_index; i_P++) {
+ vector<int> 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; ielem<nb_intersecting_elems; ielem++) {
+ //BBTree structure returns numbers between 0 and n-1
+ int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+ int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL, i_S+1)
+ %100;
+ double surf=intersector->intersectCells(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<<endl<<"Printing intersection areas:"<<endl<<endl;
+// cout<<"(source cell, target cell): intersection areas"<<endl;
+// double total=0.0;
+// double total_interm=0.0;
+// int nb_result_areas = result_areas.size();
+// for (int i=0; i< nb_result_areas; i++) {
+// map<int,double>::iterator surface;
+// total_interm=0.0;
+// for (surface=result_areas[i].begin(); surface
+// !=result_areas[i].end(); surface++) {
+// cout<<" ("<<i+1<<" , "<<(*surface).first<<")"<<" : "<<(*surface).second<<endl;
+// total_interm +=(*surface).second;
+// }
+// cout<< " elem " << i+1 << " area= " << total_interm << endl;
+// total+=total_interm;
+// }
+// cout << "total area "<<total<<endl;
+// }
+// cout<< "Filtering time= " << end_filtering-start_filtering<< endl;
+// cout<< "Intersection time= " << end_intersection-start_intersection
+// << endl;
+// long global_end =clock();
+// cout<< "Global time= " << global_end - global_start << endl;
+// }
+
+
+}
+
+//
+///** \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
+// *
+// */
+//MEDMEM_Matrix<double,1>* Interpolation3DSurf::interpolate(
+// 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"
+// <<endl;
+// string* Type_S = myMesh_S.getCellTypeNames(MED_EN::MED_CELL);
+// string* Type_P = myMesh_P.getCellTypeNames(MED_EN::MED_CELL);
+//
+// if ((Type_S[0] != "MED_TRIA3" && Type_S[0] != "MED_QUAD4") || (Type_P[0]
+// != "MED_TRIA3" && Type_P[0] != "MED_QUAD4"))
+// cout
+// <<"Interpolation3DSurf::interpolateMeshes: one of the two meshes is neither linear triangular nor linear quadratic"
+// <<endl;
+// //VB
+// delete[] Type_S;
+// delete[] Type_P;
+// int nbMaille_S =myMesh_S.getNumberOfElements(MED_EN::MED_CELL,
+// MED_EN::MED_ALL_ELEMENTS);
+// int nbMaille_P =myMesh_P.getNumberOfElements(MED_EN::MED_CELL,
+// MED_EN::MED_ALL_ELEMENTS);
+//
+// /**************************************************/
+// /* Search the characteristic size of the meshes */
+// /**************************************************/
+//
+// vector<vector<double> > BoxS=myMesh_S.getBoundingBox();
+// vector<vector<double> > 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<<endl;
+// cout<<" - Characteristic size of the target mesh: "
+// << DimCaracteristic_P<<endl;
+// cout << "Interpolation3DSurf::computation of the intersections"<<endl;
+// }
+//
+// PlanarIntersector* intersector;
+//
+// switch (_Intersection_type) {
+// case MEDMEM::Triangulation:
+// intersector=new TriangulationIntersector<3>(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<double> 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<int,double> MAP_init;
+//
+// MEDMEM_Matrix<double,1>* matrix = new MEDMEM_Matrix<double,1> (nbMaille_P);
+// //vector<map<int,double> > 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<NumberOfCellsTypes_P; itype++) {
+// int nb_nodesP=types[itype]%100;
+// int nbelem_type = myMesh_P.getNumberOfElements(MED_EN::MED_CELL,
+// types[itype]);
+// type_max_index += nbelem_type;
+// for (i_P=type_min_index; i_P<type_max_index; i_P++) {
+// vector<int> 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; ielem<nb_intersecting_elems; ielem++) {
+// //BBTree structure returns numbers between 0 and n-1
+// int i_S=intersecting_elems[ielem]; //MN: Global number of cell ?
+// int nb_nodesS=myMesh_S.getElementType(MED_EN::MED_CELL, i_S+1)
+// %100;
+// double surf=intersector->intersectCells(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;
+//}
+}
+;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /* 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,
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /* 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 */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
{
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 ;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /*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. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
// 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<taille/2;i++)
+ for(uint i=0;i<taille/2;i++)
{
if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))<absolute_precision)
isPresent=true;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/* fonction qui rajoute les sommet du triangle P dans le vecteur V */
- /* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans */
+ /* si ceux-ci sont compris dans le triangle S et ne sont pas d�j� dans */
/* V. */
/*sommets de P: P_1, P_2, P_3 */
/*sommets de S: P_4, P_5, P_6 */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */
/* . Si l'intersection est non nulle et si celle-ci n'est */
- /* n'est pas déjà contenue dans Vect on la rajoute à Vect */
+ /* n'est pas d�j� contenue dans Vect on la rajoute � Vect */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
inline void inters_de_segment(const double * P_1,const double * P_2,
{
- // calcul du déterminant de P_1P_2 et P_3P_4.
+ // calcul du d�terminant de P_1P_2 et P_3P_4.
double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]);
double absolute_precision = dim_caracteristic*precision;
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/* calcul l'intersection de deux triangles */
/* P_1, P_2, P_3: sommets du premier triangle */
- /* P_4, P_5, P_6: sommets du deuxième triangle */
+ /* P_4, P_5, P_6: sommets du deuxi�me triangle */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
int N_2=MaS._Connect[3*(i_S-1)+1];
int N_3=MaS._Connect[3*(i_S-1)+2];
vector<int> 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;
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)
{
}
}
- for(int i=0;i<V_23.size();i++)
+ for(uint i=0;i<V_23.size();i++)
{
if(V_23[i]!=i_S)
{
}
}
- for(int i=0;i<V_13.size();i++)
+ for(uint i=0;i<V_13.size();i++)
{
if(V_13[i]!=i_S)
{
}
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
- /* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer */
+ /* fonction pour v�rifier qu'un n�de maille n'a pas d�ja �t� consid�rer */
/* dans le vecteur et le rajouter au vecteur sinon. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /* fonction pour reconstituer un polygone convexe à partir */
+ /* fonction pour reconstituer un polygone convexe � partir */
/* d'un nuage de point. */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
- /* fonction qui écrit les résultats d'annelise dans un fichier: */
+ /* fonction qui �crit les r�sultats d'annelise dans un fichier: */
/* pour chaque maille du maillage 1 on stoque la fraction volumique */
- /* considéré lors de l'algorithme */
+ /* consid�r� lors de l'algorithme */
/*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
inline void WriteFractSurfInFile(MEDMEM::FIELD<double>& 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");
// 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 {
}
// 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];
}
{
// count negative as face has reversed orientation
totalVolume -= _volumes[key2];
- }
+ }
}
+
break;
default:
ConvexIntersector.hxx\
ConvexIntersector.cxx\
GenericIntersector.hxx\
-GenericIntersector.cxx
+GenericIntersector.cxx\
+Remapper.hxx
# Libraries targets
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 =
--- /dev/null
+#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<double,1>;
+ 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<double>& field_source,MEDMEM::FIELD<double>& 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<double>* target_volumes = getSupportVolumes(*field_target.getSupport());
+ int nbelem_target=field_target.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
+
+ double* value_target = const_cast<double*> (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<double>* 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");
+ }
+}
+
+
+}
--- /dev/null
+#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<double>& field_source, MEDMEM::FIELD<double>& field_target);
+ void setOptionDouble(const std::string& key, double value);
+ void setOptionInt(const std::string& key, int value);
+private :
+ MEDMEM::MEDMEM_Matrix<double,1>* _matrix;
+ MEDMEM::FIELD<double>* getSupportVolumes(const MEDMEM::SUPPORT& support);
+
+};
+
+}
+
+#endif /*REMAPPER_HXX_*/
Interpolation3DTestSuite.hxx \
SingleElementTetraTests.hxx \
MultiElementTetraTests.hxx \
- HexaTests.hxx\
- MeshTestToolkit.hxx
+ HexaTests.hxx \
+ MeshTestToolkit.hxx \
+ BBTreeTest.hxx
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 =
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.
#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.
#include "MeshTestToolkit.hxx"
#include "MEDMEM_Mesh.hxx"
+#include "Interpolation3D.txx"
#include <iostream>
#include <map>
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)
#include "Interpolation3D.hxx"
+#include "Interpolation3D.txx"
#include "MeshTestToolkit.hxx"
#include "Log.hxx"
#include "VectorUtils.hxx"
LOG(1, "Target mesh has " << numTargetElems << " elements");
- m = _interpolator->interpolateMeshes(sMesh, tMesh);
+ _interpolator->interpolateMeshes(sMesh, tMesh,m);
std::pair<int, int> eff = countNumberOfMatrixEntries(m);
LOG(1, eff.first << " of " << numTargetElems * numSrcElems << " intersections calculated : ratio = "
#include "MultiElementTetraTests.hxx"
#include "SingleElementTetraTests.hxx"
#include "HexaTests.hxx"
+#include "BBTreeTest.hxx"
using namespace INTERP_TEST;
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"
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;
