template<class MyMeshType, class MatrixType>
int toIntegralUniform(const MyMeshType& meshS, MatrixType& result, const char *method) { return fromToIntegralUniform(true,meshS,result,method); }
static void checkAndSplitInterpolationMethod(const char *method, std::string& srcMeth, std::string& trgMeth) throw(INTERP_KERNEL::Exception);
+ template<class MyMeshType>
+ static double CalculateCharacteristicSizeOfMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, const int printLevel);
protected:
template<class MyMeshType, class MatrixType>
int fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method);
#include "Interpolation.hxx"
#include "IntegralUniformIntersector.hxx"
#include "IntegralUniformIntersector.txx"
+#include "VectorUtils.hxx"
namespace INTERP_KERNEL
{
srcMeth=methodC.substr(0,2);
trgMeth=methodC.substr(2);
}
+
+ template<class TrueMainInterpolator>
+ template<class MyMeshType>
+ double Interpolation<TrueMainInterpolator>::CalculateCharacteristicSizeOfMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, const int printLevel)
+ {
+ static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
+
+ long nbMailleS=myMeshS.getNumberOfElements();
+ long nbMailleT=myMeshT.getNumberOfElements();
+
+ /**************************************************/
+ /* Search the characteristic size of the meshes */
+ /**************************************************/
+
+ double BoxS[2*SPACEDIM]; myMeshS.getBoundingBox(BoxS);
+ double BoxT[2*SPACEDIM]; myMeshT.getBoundingBox(BoxT);
+ double diagonalS,dimCaracteristicS=std::numeric_limits<double>::max();
+ if(nbMailleS!=0)
+ {
+ diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
+ dimCaracteristicS=diagonalS/nbMailleS;
+ }
+ double diagonalT,dimCaracteristicT=std::numeric_limits<double>::max();
+ if(nbMailleT!=0)
+ {
+ diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
+ dimCaracteristicT=diagonalT/nbMailleT;
+ }
+ if (printLevel>=1)
+ {
+ std::cout << " - Characteristic size of the source mesh : " << dimCaracteristicS << std::endl;
+ std::cout << " - Characteristic size of the target mesh: " << dimCaracteristicT << std::endl;
+ }
+
+ return std::min(dimCaracteristicS, dimCaracteristicT);
+
+ }
+
}
#endif
#ifndef __INTERPOLATION3D2D_HXX__
#define __INTERPOLATION3D2D_HXX__
+#include <set>
+#include <map>
+
#include "Interpolation.hxx"
#include "NormalizedUnstructuredMesh.hxx"
#include "InterpolationOptions.hxx"
class Interpolation3D2D : public Interpolation<Interpolation3D2D>
{
public:
+ typedef typename std::map<int,std::set<int> > DuplicateFacesType;
+
Interpolation3D2D();
Interpolation3D2D(const InterpolationOptions& io);
- template<class MyMeshType, class MatrixType>
- int interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method);
+ template<class MyMeshType, class MyMatrixType>
+ int interpolateMeshes(const MyMeshType& srcMesh,
+ const MyMeshType& targetMesh,
+ MyMatrixType& matrix,
+ const char *method);
+ DuplicateFacesType retrieveDuplicateFaces() const
+ {
+ return _duplicate_faces;
+ }
private:
SplittingPolicy _splitting_policy;
+ DuplicateFacesType _duplicate_faces;
};
}
#include "PolyhedronIntersectorP1P1.txx"
#include "PointLocator3DIntersectorP1P1.txx"
#include "Log.hxx"
-/// 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.txx"
-#include "RegionNode.hxx"
-#include <stack>
-
-#else // use BBTree class
#include "BBTree.txx"
-#endif
-
namespace INTERP_KERNEL
{
/**
* @param srcMesh 3-dimensional source mesh
* @param targetMesh 3-dimesional target mesh, containing only tetraedra
- * @param result matrix in which the result is stored
+ * @param matrix matrix in which the result is stored
*
*/
- template<class MyMeshType, class MatrixType>
- int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method)
+ template<class MyMeshType, class MyMatrixType>
+ int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh,
+ const MyMeshType& targetMesh,
+ MyMatrixType& matrix,
+ const char *method)
{
typedef typename MyMeshType::MyConnType ConnType;
// create MeshElement objects corresponding to each element of the two meshes
std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
std::map<MeshElement<ConnType>*, int> indices;
+ DuplicateFacesType intersectFaces;
for(unsigned long i = 0 ; i < numSrcElems ; ++i)
srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
for(unsigned long i = 0 ; i < numTargetElems ; ++i)
targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
- Intersector3D<MyMeshType,MatrixType>* intersector=0;
+ Intersector3D<MyMeshType,MyMatrixType>* intersector=0;
std::string methC = InterpolationOptions::filterInterpolationMethod(method);
+ const double dimCaracteristic = CalculateCharacteristicSizeOfMeshes(srcMesh, targetMesh, InterpolationOptions::getPrintLevel());
if(methC=="P0P0")
{
switch(InterpolationOptions::getIntersectionType())
{
case Triangulation:
- intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MatrixType>(targetMesh, srcMesh,
- getSplittingPolicy());
+ intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>(targetMesh,
+ srcMesh,
+ dimCaracteristic,
+ getPrecision(),
+ intersectFaces,
+ getSplittingPolicy());
break;
default:
throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangulation.");
else
throw Exception("Invalid method choosed must be in \"P0P0\".");
// create empty maps for all source elements
- result.resize(intersector->getNumberOfRowsOfResMatrix());
-
- // TODO DP : #ifdef USE_RECURSIVE_BBOX_FILTER : voir Interpolation3D2D.txx
+ matrix.resize(intersector->getNumberOfRowsOfResMatrix());
// create BBTree structure
// - get bounding boxes
srcElemIdx[i] = srcElems[i]->getIndex();
}
-#if 0//dp
- BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems);
-#else
BBTree<3,ConnType> tree(bboxes, srcElemIdx, 0, numSrcElems, 0.);
-#endif
// for each target element, get source elements with which to calculate intersection
// - calculate intersection by calling intersectCells
tree.getIntersectingElems(targetBox, intersectElems);
if ( !intersectElems.empty() )
- intersector->intersectCells(targetIdx,intersectElems,result);
+ intersector->intersectCells(targetIdx, intersectElems, matrix);
+
}
delete [] bboxes;
delete [] srcElemIdx;
+ DuplicateFacesType::iterator iter;
+ for (iter = intersectFaces.begin(); iter != intersectFaces.end(); ++iter)
+ {
+ if (iter->second.size() > 1)
+ {
+ _duplicate_faces.insert(std::make_pair(iter->first, iter->second));
+ }
+ }
+
// free allocated memory
int ret=intersector->getNumberOfColsOfResMatrix();
// Main function to interpolate triangular and quadratic meshes
template<class MyMeshType, class MatrixType>
int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method);
+
public:
bool doRotate() const { return asLeafInterpPlanar().doRotate(); }
double medianPlane() const { return asLeafInterpPlanar().medianPlane(); }
InterpolationOptions::setIntersectionType(intersectionType);
InterpolationOptions::setOrientation(orientation);
}
-
-
+
/** \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
/***********************************************************/
long nbMailleS=myMeshS.getNumberOfElements();
- long nbMailleT=myMeshT.getNumberOfElements();
/**************************************************/
/* Search the characteristic size of the meshes */
/**************************************************/
- double BoxS[2*SPACEDIM]; myMeshS.getBoundingBox(BoxS);
- double BoxT[2*SPACEDIM]; myMeshT.getBoundingBox(BoxT);
- double diagonalS,dimCaracteristicS=std::numeric_limits<double>::max();
- if(nbMailleS!=0)
- {
- diagonalS=getDistanceBtw2Pts<SPACEDIM>(BoxS+SPACEDIM,BoxS);
- dimCaracteristicS=diagonalS/nbMailleS;
- }
- double diagonalT,dimCaracteristicT=std::numeric_limits<double>::max();
- if(nbMailleT!=0)
- {
- diagonalT=getDistanceBtw2Pts<SPACEDIM>(BoxT+SPACEDIM,BoxT);
- dimCaracteristicT=diagonalT/nbMailleT;
- }
-
- _dim_caracteristic=std::min(dimCaracteristicS, dimCaracteristicT);
- if (InterpolationOptions::getPrintLevel()>=1)
+ int printLevel = InterpolationOptions::getPrintLevel();
+ _dim_caracteristic = CalculateCharacteristicSizeOfMeshes(myMeshS, myMeshT, printLevel);
+ if (printLevel>=1)
{
- std::cout << " - Characteristic size of the source mesh : " << dimCaracteristicS << std::endl;
- std::cout << " - Characteristic size of the target mesh: " << dimCaracteristicT << std::endl;
std::cout << "InterpolationPlanar::computation of the intersections" << std::endl;
}
#include "INTERPKERNELDefines.hxx"
#include "InterpKernelException.hxx"
-#if 1//dp
-#include "VectorUtils.hxx"
-#endif
#include "NormalizedUnstructuredMesh.hxx"
* the source elements.
*
*/
- template<class MyMeshType, class MyMatrix>
- class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0<MyMeshType,MyMatrix>
- {
+ template<class MyMeshType, class MyMatrixType>
+ class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0<MyMeshType,MyMatrixType>
+ {
+ typedef typename std::map<int,std::set<int> > DuplicateFacesType;
+
public:
static const int SPACEDIM=MyMeshType::MY_SPACEDIM;
static const int MESHDIM=MyMeshType::MY_MESHDIM;
typedef typename MyMeshType::MyConnType ConnType;
static const NumberingPolicy numPol=MyMeshType::My_numPol;
+
public:
- Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, const MyMeshType& srcMesh,
+ Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh,
+ const MyMeshType& srcMesh,
+ const double dimCaracteristic,
+ const double precision,
+ DuplicateFacesType& intersectFaces,
SplittingPolicy policy = PLANAR_FACE_5);
~Polyhedron3D2DIntersectorP0P0();
- void intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res);
+ void intersectCells(ConnType targetCell,
+ const std::vector<ConnType>& srcCells,
+ MyMatrixType& matrix);
private:
void releaseArrays();
SplitterTetra2<MyMeshType> _split;
+ double _dim_caracteristic;
+ double _precision;
+
+ DuplicateFacesType& _intersect_faces;
+
};
}
* @param srcMesh mesh containing the source elements
* @param policy splitting policy to be used
*/
- template<class MyMeshType, class MyMatrix>
- Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, const MyMeshType& srcMesh,
- SplittingPolicy policy)
- : Intersector3DP0P0<MyMeshType,MyMatrix>(targetMesh,srcMesh),_split(targetMesh,srcMesh,policy)
+ template<class MyMeshType, class MyMatrixType>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh,
+ const MyMeshType& srcMesh,
+ const double dimCaracteristic,
+ const double precision,
+ DuplicateFacesType& intersectFaces,
+ SplittingPolicy policy)
+ : Intersector3DP0P0<MyMeshType,MyMatrixType>(targetMesh,srcMesh),
+ _split(targetMesh,srcMesh,policy),
+ _dim_caracteristic(dimCaracteristic),
+ _precision(precision),
+ _intersect_faces(intersectFaces)
{
}
* Liberates the SplitterTetra objects and potential sub-node points that have been allocated.
*
*/
- template<class MyMeshType, class MyMatrix>
- Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::~Polyhedron3D2DIntersectorP0P0()
+ template<class MyMeshType, class MyMatrixType>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::~Polyhedron3D2DIntersectorP0P0()
{
releaseArrays();
}
- template<class MyMeshType, class MyMatrix>
- void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::releaseArrays()
+ template<class MyMeshType, class MyMatrixType>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::releaseArrays()
{
for(typename std::vector< SplitterTetra<MyMeshType>* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
delete *iter;
* @param srcCells in C mode.
*
*/
- template<class MyMeshType, class MyMatrix>
- void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
+ template<class MyMeshType, class MyMatrixType>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrixType>::intersectCells(ConnType targetCell,
+ const std::vector<ConnType>& srcCells,
+ MyMatrixType& matrix)
{
- int nbOfNodesT=Intersector3D<MyMeshType,MyMatrix>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell));
+ int nbOfNodesT=Intersector3D<MyMeshType,MyMatrixType>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell));
releaseArrays();
_split.splitTargetCell(targetCell,nbOfNodesT,_tetra);
for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
{
double surface = 0.;
- std::set<TriangleFaceKey> listOfTetraFacesTreated;
+ std::multiset<TriangleFaceKey> listOfTetraFacesTreated;
+ std::set<TriangleFaceKey> listOfTetraFacesColinear;
// calculate the coordinates of the nodes
const NumberingPolicy numPol=MyMeshType::My_numPol;
typename MyMeshType::MyConnType cellSrc = *iterCellS;
int cellSrcIdx = OTT<ConnType,numPol>::indFC(cellSrc);
- NormalizedCellType normCellType=Intersector3D<MyMeshType,MyMatrix>::_src_mesh.getTypeOfElement(cellSrcIdx);
+ NormalizedCellType normCellType=Intersector3D<MyMeshType,MyMatrixType>::_src_mesh.getTypeOfElement(cellSrcIdx);
const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
- const MyMeshType& _src_mesh = Intersector3D<MyMeshType,MyMatrix>::_src_mesh;
+ const MyMeshType& _src_mesh = Intersector3D<MyMeshType,MyMatrixType>::_src_mesh;
unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(cellSrcIdx) : cellModelCell.getNumberOfNodes();
int *polyNodes=new int[nbOfNodes4Type];
double **polyCoords = new double*[nbOfNodes4Type];
}
for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
- surface += (*iter)->intersectSourceFace(normCellType, nbOfNodes4Type, polyNodes, polyCoords, listOfTetraFacesTreated);
- if(surface!=0.)
- res[targetCell].insert(std::make_pair(OTT<ConnType,numPol>::indFC(*iterCellS), surface));
- listOfTetraFacesTreated.clear();
+ surface += (*iter)->intersectSourceFace(normCellType,
+ nbOfNodes4Type,
+ polyNodes,
+ polyCoords,
+ _dim_caracteristic,
+ _precision,
+ listOfTetraFacesTreated,
+ listOfTetraFacesColinear);
+
+ if(surface!=0.) {
+
+ matrix[targetCell].insert(std::make_pair(cellSrcIdx, surface));
+
+ bool isSrcFaceColinearWithFaceOfTetraTargetCell = false;
+ std::set<TriangleFaceKey>::iterator iter;
+ for (iter = listOfTetraFacesColinear.begin(); iter != listOfTetraFacesColinear.end(); ++iter)
+ {
+ if (listOfTetraFacesTreated.count(*iter) != 1)
+ {
+ isSrcFaceColinearWithFaceOfTetraTargetCell = false;
+ break;
+ }
+ else
+ {
+ isSrcFaceColinearWithFaceOfTetraTargetCell = true;
+ }
+ }
+
+ if (isSrcFaceColinearWithFaceOfTetraTargetCell)
+ {
+ DuplicateFacesType::iterator intersectFacesIter = _intersect_faces.find(cellSrcIdx);
+ if (intersectFacesIter != _intersect_faces.end())
+ {
+ intersectFacesIter->second.insert(targetCell);
+ }
+ else
+ {
+ std::set<int> targetCellSet;
+ targetCellSet.insert(targetCell);
+ _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet));
+ }
+
+ }
+
+ }
+
delete[] polyNodes;
delete[] polyCoords;
#include <vector>
#include <functional>
#include <map>
+#include <set>
+//dp#include <algorithm>
namespace INTERP_KERNEL
{
const int polyNodesNbr,
const int *const polyNodes,
const double *const *const polyCoords,
- std::set<TriangleFaceKey>& listOfTetraFacesTreated);
+ const double dimCaracteristic,
+ const double precision,
+ std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
+ std::set<TriangleFaceKey>& listOfTetraFacesColinear);
double intersectTetra(const double** tetraCorners);
// member functions
inline void createAffineTransform(const double** corners);
inline void checkIsOutside(const double* pt, bool* isOutside) const;
- inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const; //dp à supprimer
+ inline void checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol = DEFAULT_ABS_TOL) const;
inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum);
+ inline void calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node);
inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key);
inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key);
inline bool isFacesCoplanar(const double *const plane_normal, const double plane_constant,
const double *const *const coordsFace);
inline double calculateIntersectionSurfaceOfCoplanarTriangles(const double *const normal,
+ const double plane_constant,
const double *const P_1, const double *const P_2, const double *const P_3,
const double *const P_4, const double *const P_5, const double *const P_6,
double dim_caracteristic, double precision);
isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
}
-#if 1//dp
- //dp à supprimer
template<class MyMeshType>
inline void SplitterTetra<MyMeshType>::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const
{
isStrictlyOutside[6] = isStrictlyOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol);
isStrictlyOutside[7] = isStrictlyOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol));
}
-#endif
/**
* Calculates the transformed node with a given global node number.
_nodes[globalNodeNum] = transformedNode;
}
+ //dp TODO DP : supprimer la précédente et commentariser
+ template<class MyMeshType>
+ inline void SplitterTetra<MyMeshType>::calculateNode2(typename MyMeshType::MyConnType globalNodeNum, const double* node)
+ {
+ double* transformedNode = new double[MyMeshType::MY_SPACEDIM];
+ assert(transformedNode != 0);
+ _t->apply(transformedNode, node);
+ _nodes[globalNodeNum] = transformedNode;
+ }
+
/**
* Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in .
* Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation.
faceNodes=new int[faceModel.getNumberOfNodes()];
cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes);
}
-
// intersect a son with the unit tetra
switch(faceType)
{
return std::fabs(1.0 / _t->determinant() * totalVolume) ;
}
-#if 1//dp TODO DP : à commenter
+ //dp TODO DP : à commenter
template<class MyMeshType>
double SplitterTetra<MyMeshType>::calculateIntersectionSurfaceOfCoplanarTriangles(const double *const plane_normal,
+ const double plane_constant,
const double *const P_1, const double *const P_2, const double *const P_3,
const double *const P_4, const double *const P_5, const double *const P_6,
double dim_caracteristic, double precision)
inter2,
dim_caracteristic, precision);
ConnType nb_inter=((ConnType)inter2.size())/2;
- if(nb_inter >3) inter2=reconstruct_polygon(inter2);
double surface = 0.;
- double normal[3];
- for(ConnType i = 1; i<nb_inter-1; i++)
- {
- INTERP_KERNEL::crossprod<2>(&inter2[0],&inter2[2*i],&inter2[2*(i+1)],normal);
- surface +=0.5*fabs(normal[0]);
- }
- return surface;
-#if 0//dp
- int nb_inter = inter2.size();
+ if(nb_inter >3) inter2=reconstruct_polygon(inter2);
if (nb_inter > 0)
{
- inter3.resize(3 * nb_inter / 2);
+ std::vector<double> inter3;
+ inter3.resize(3 * nb_inter);
// Map 2D intersections back to the 3D triangle space.
if (maxNormal == 0)
{
- double invNX = ((double) 1.) / normal[0];
+ double invNX = ((double) 1.) / plane_normal[0];
for (i = 0; i < nb_inter; i++)
{
inter3[3 * i + 1] = inter2[2 * i];
inter3[3 * i + 2] = inter2[2 * i + 1];
- inter3[3 * i] = invNX * (constant - normal[1] * inter3[3 * i + 1] - normal[2] * inter3[3 * i + 2]);
+ inter3[3 * i] = invNX * (plane_constant - plane_normal[1] * inter3[3 * i + 1] - plane_normal[2] * inter3[3 * i + 2]);
}
}
else if (maxNormal == 1)
{
- double invNY = ((double) 1.) / normal[1];
+ double invNY = ((double) 1.) / plane_normal[1];
for (i = 0; i < nb_inter; i++)
{
inter3[3 * i] = inter2[2 * i];
inter3[3 * i + 2] = inter2[2 * i + 1];
- inter3[3 * i + 1] = invNY * (constant - normal[0] * inter3[3 * i] - normal[2] * inter3[3 * i + 2]);
+ inter3[3 * i + 1] = invNY * (plane_constant - plane_normal[0] * inter3[3 * i] - plane_normal[2] * inter3[3 * i + 2]);
}
}
else
{
- double invNZ = ((double) 1.) / normal[2];
+ double invNZ = ((double) 1.) / plane_normal[2];
for (i = 0; i < nb_inter; i++)
{
inter3[3 * i] = inter2[2 * i];
inter3[3 * i + 1] = inter2[2 * i + 1];
- inter3[3 * i + 2] = invNZ * (constant - normal[0] * inter3[3 * i] - normal[1] * inter3[3 * i + 1]);
+ inter3[3 * i + 2] = invNZ * (plane_constant - plane_normal[0] * inter3[3 * i] - plane_normal[1] * inter3[3 * i + 1]);
}
}
+ surface = polygon_area<3>(inter3);
}
-#endif
+ return surface;
}
template<class MyMeshType>
- bool SplitterTetra<MyMeshType>::isFacesCoplanar(const double *const plane_normal, const double plane_constant,
+ bool SplitterTetra<MyMeshType>::isFacesCoplanar(const double *const plane_normal,
+ const double plane_constant,
const double *const *const coordsFace)
{
// Compute the signed distances of triangle vertices to the plane. Use an epsilon-thick plane test.
return false;
}
-#endif
-
-
// TODO DP : adapter les commentaires suivants. _volume => _surface ?
/**
* Calculates the volume of intersection of an element in the source mesh and the target element.
const int polyNodesNbr,
const int *const polyNodes,
const double *const *const polyCoords,
- std::set<TriangleFaceKey>& listOfTetraFacesTreated)
+ const double dimCaracteristic,
+ const double precision,
+ std::multiset<TriangleFaceKey>& listOfTetraFacesTreated,
+ std::set<TriangleFaceKey>& listOfTetraFacesColinear)
{
typedef typename MyMeshType::MyConnType ConnType;
- double totalVolume = 0.0;
-
+ double totalSurface = 0.0;
-#if 0
- // Is src element coplanar with one of the tetra faces ?
-
- double srcNormal[3];
-#if 0//dp
- const double* points[3];
- for(int i = 0 ; i < 3 ; ++i)
- {
- points[i] = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*cellNodes[i];
- }
- calculateNormalForTria(points[0],points[1],points[2], srcNormal);
-#else
- calculateNormalForPolyg((const double **)coordsPoly, nbOfNodes4Type, srcNormal);
-#endif
-
- double faceNormal[3];
- double crossNormals[3];
- const int tetraFacesNodesConn[4][3] = {{0,1,2},{0,2,3},{0,3,1},{1,2,3}};
- for (int iTetraFace = 0; iTetraFace < 4; ++iTetraFace)
- {
- const int *const tetraFaceNodesConn = tetraFacesNodesConn[iTetraFace];
- const double *const coordsTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM;
- const double *const coordsTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM;
- const double *const coordsTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM;
- calculateNormalForTria(coordsTriNode1, coordsTriNode2, coordsTriNode3, faceNormal);
- cross(srcNormal, faceNormal, crossNormals);
- if (epsilonEqual(norm(crossNormals), 0.))
- {
- // Les faces sont sur des plans parallèles
- double area[3];
- int nbTria = nbOfNodes4Type - 2; // split polygon into nbTria triangles
- for (int iTri = 0; iTri < nbTria; ++iTri)
- {
- std::vector<double> inter;
- INTERP_KERNEL::intersec_de_triangle(coordsPoly[0], coordsPoly[1 + iTri], coordsPoly[2 + iTri],
- coordsTriNode1, coordsTriNode2, coordsTriNode3,
- inter,
- 1., //dp inter, PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
- DEFAULT_ABS_TOL); //dp PlanarIntersector<MyMeshType,MyMatrix>::_precision);
- ConnType nb_inter=((ConnType)inter.size())/2;
- if(nb_inter >3) inter=reconstruct_polygon(inter);
- for(ConnType i = 1; i<nb_inter-1; i++)
- {
- INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
- totalVolume +=0.5*fabs(area[0]);
- }
- }
- }
- }
-#endif
-
- //{ could be done on outside?
// check if we have planar tetra element
if(_t->determinant() == 0.0)
{
bool isTargetOutside = false;
// calculate the coordinates of the nodes
-#if 0//dp
- int *polyNodes=new int[polyNodesNbr];
-#endif
for(int i = 0;i<(int)polyNodesNbr;++i)
{
-#if 0//dp
- // we could store mapping local -> global numbers too, but not sure it is worth it
- const int globalNodeNum = getGlobalNumberOfNode(i, OTT<ConnType,numPol>::indFC(element), _src_mesh);
- polyNodes[i]=globalNodeNum;
-#else
const int globalNodeNum = polyNodes[i];
-#endif
if(_nodes.find(globalNodeNum) == _nodes.end())
{
//for(HashMap< int , double* >::iterator iter3=_nodes.begin();iter3!=_nodes.end();iter3++)
// std::cout << (*iter3).first << " ";
//std::cout << std::endl << "*** " << globalNodeNum << std::endl;
- calculateNode(globalNodeNum);
+ calculateNode2(globalNodeNum, polyCoords[i]);
}
checkIsStrictlyOutside(_nodes[globalNodeNum], isStrictlyOutside);
_conn[tetraFaceNodesConn[2]]);
if (listOfTetraFacesTreated.find(key) == listOfTetraFacesTreated.end())
{
- listOfTetraFacesTreated.insert(key);
const double * const coordsTetraTriNode1 = _coords + tetraFaceNodesConn[0] * MyMeshType::MY_SPACEDIM;
const double * const coordsTetraTriNode2 = _coords + tetraFaceNodesConn[1] * MyMeshType::MY_SPACEDIM;
const double * const coordsTetraTriNode3 = _coords + tetraFaceNodesConn[2] * MyMeshType::MY_SPACEDIM;
int nbrPolyTri = polyNodesNbr - 2; // split polygon into nbrPolyTri triangles
for (int iTri = 0; iTri < nbrPolyTri; ++iTri)
{
- totalVolume += calculateIntersectionSurfaceOfCoplanarTriangles(plane_normal,
- polyCoords[0],
- polyCoords[1 + iTri],
- polyCoords[2 + iTri],
- coordsTetraTriNode1,
- coordsTetraTriNode2,
- coordsTetraTriNode3,
- 1., //dp PlanarIntersector<MyMeshType,MyMatrix>::_dim_caracteristic,
- DEFAULT_ABS_TOL); //dp PlanarIntersector<MyMeshType,MyMatrix>::_precision);
+ double volume = calculateIntersectionSurfaceOfCoplanarTriangles(plane_normal,
+ plane_constant,
+ polyCoords[0],
+ polyCoords[1 + iTri],
+ polyCoords[2 + iTri],
+ coordsTetraTriNode1,
+ coordsTetraTriNode2,
+ coordsTetraTriNode3,
+ dimCaracteristic,
+ precision);
+ if (!epsilonEqual(volume, 0.))
+ {
+ totalSurface += volume;
+ listOfTetraFacesColinear.insert(key);
+ }
}
}
}
+ listOfTetraFacesTreated.insert(key);
}
}
else
{
TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
calculateSurface(tri, key);
- totalVolume += _volumes[key];
+ totalSurface += _volumes[key];
}
else
{
// count negative as face has reversed orientation
- totalVolume -= _volumes[key];
+ totalSurface -= _volumes[key];
}
}
break;
- case NORM_QUAD4: //dp pas d'intérêt vis-à-vis du case suivant (à regrouper les deux cases)
+ case NORM_QUAD4:
// simple triangulation of faces along a diagonal :
//
{
TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1]], _nodes[polyNodes[2]]);
calculateSurface(tri, key1);
- totalVolume += _volumes[key1];
+ totalSurface += _volumes[key1];
}
else
{
// count negative as face has reversed orientation
- totalVolume -= _volumes[key1];
+ totalSurface -= _volumes[key1];
}
// local nodes 1, 3, 4
{
TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[2]], _nodes[polyNodes[3]]);
calculateSurface(tri, key2);
- totalVolume += _volumes[key2];
+ totalSurface += _volumes[key2];
}
else
{
// count negative as face has reversed orientation
- totalVolume -= _volumes[key2];
+ totalSurface -= _volumes[key2];
}
}
break;
TransformedTriangle tri(_nodes[polyNodes[0]], _nodes[polyNodes[1 + iTri]],
_nodes[polyNodes[2 + iTri]]);
calculateSurface(tri, key);
- totalVolume += _volumes[key];
+ totalSurface += _volumes[key];
}
else
{
- totalVolume -= _volumes[key];
+ totalSurface -= _volumes[key];
}
}
}
// reset if it is very small to keep the matrix sparse
// is this a good idea?
- if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT))
+ if(epsilonEqual(totalSurface, 0.0, SPARSE_TRUNCATION_LIMIT))
{
- totalVolume = 0.0;
+ totalSurface = 0.0;
}
- LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant());
+ LOG(2, "Volume = " << totalSurface << ", det= " << _t->determinant());
- return totalVolume;
+ return totalSurface;
}
/**
int conn[4];
for(int j = 0; j < 4; ++j)
{
+#if 0
nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_5[4*i+j] ],conn[j]);
+#else
+ conn[j] = subZone[ SPLIT_NODES_5[4*i+j] ];
+ nodes[j] = getCoordsOfSubNode(conn[j]);
+#endif
}
SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
tetra.push_back(t);
int conn[4];
for(int j = 0; j < 4; ++j)
{
-#if 1//dp
conn[j] = subZone[SPLIT_NODES_6[4*i+j]];
nodes[j] = getCoordsOfSubNode(conn[j]);
-#else
- nodes[j] = getCoordsOfSubNode2(subZone[SPLIT_NODES_6[4*i+j]], conn[j]);
-#endif
}
SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
tetra.push_back(t);
const double* nodes[4];
int conn[4];
// get the cell center
- nodes[0] = getCoordsOfSubNode2(14,conn[0]);
-
+ conn[0] = 14;
+ nodes[0] = getCoordsOfSubNode(conn[0]);
+
for(int faceCenterNode = 8; faceCenterNode < 14; ++faceCenterNode)
{
// get the face center
- nodes[1] = getCoordsOfSubNode2(faceCenterNode,conn[1]);
+ conn[1] = faceCenterNode;
+ nodes[1] = getCoordsOfSubNode(conn[1]);
for(int j = 0; j < 4; ++j)
{
const int row = 4*(faceCenterNode - 8) + j;
- nodes[2] = getCoordsOfSubNode2(TETRA_EDGES[2*row],conn[2]);
- nodes[3] = getCoordsOfSubNode2(TETRA_EDGES[2*row + 1],conn[3]);
-
+ conn[2] = TETRA_EDGES[2*row];
+ conn[3] = TETRA_EDGES[2*row + 1];
+ nodes[2] = getCoordsOfSubNode(conn[2]);
+ nodes[3] = getCoordsOfSubNode(conn[3]);
+
SplitterTetra<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes, conn);
tetra.push_back(t);
}
#include "TransformedTriangle.hxx"
#include "VectorUtils.hxx"
-#if 1//dp
#include "TetraAffineTransform.hxx"
-#endif
+//dp #include "InterpolationUtils.hxx"
#include <iostream>
#include <fstream>
#include <cassert>
for(int i = 0 ; i < nbPoints ; ++i)
{
const double *const ptCurr = _polygonA[i]; // pt "i"
- const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt m == pt 0)
+ const double *const ptNext = _polygonA[(i + 1) % nbPoints]; // pt "i+1" (pt nbPoints == pt 0)
cross(ptCurr, ptNext, pdt);
add(pdt, sum);
};
double sign = uv_xy[0] * uv_xy[3] - uv_xy[1] * uv_xy[2];
+ if(epsilonEqual(sign, 0.))
+ {
+ sign = 0.;
+ }
return (sign < 0.) ? -1 : (sign > 0.) ? 1 : 0;
}
namespace INTERP_KERNEL
{
-#if 1//dp
class TetraAffineTransform;
-#endif
/** \class TransformedTriangle
* \brief Class representing one of the faces of the triangulated source polyhedron after having been transformed
return ss.str();
}
-
- //dp : à supprimer
/**
* Subtracts two a double[3] - vectors and store the result in res
*
namespace MEDMEM {
class MESH;
-};
+}
namespace INTERP_TEST
{
#include "TransformedTriangleIntersectTest.hxx"
#include "TransformedTriangleTest.hxx"
#include "UnitTetraIntersectionBaryTest.hxx"
+#include "UnitTetra3D2DIntersectionTest.hxx"
#ifdef DISABLE_MICROMED
#include "HexaTests.hxx"
CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleIntersectTest );
CPPUNIT_TEST_SUITE_REGISTRATION( TransformedTriangleTest );
CPPUNIT_TEST_SUITE_REGISTRATION( UnitTetraIntersectionBaryTest );
+CPPUNIT_TEST_SUITE_REGISTRATION( UnitTetra3D2DIntersectionTest );
#ifdef DISABLE_MICROMED
CPPUNIT_TEST_SUITE_REGISTRATION( HexaTests );
{
CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest);
//MEDCouplingBasicsTest1.cxx
-#if 0//dp
CPPUNIT_TEST( testArray );
CPPUNIT_TEST( testArray2 );
CPPUNIT_TEST( testArray3 );
CPPUNIT_TEST( test2DInterpP1P0PL_2 );
CPPUNIT_TEST( test2DInterpP1P1_1 );
CPPUNIT_TEST( test2DInterpP1P1PL_1 );
-#endif
CPPUNIT_TEST( test3DSurfInterpP0P0_1 );
-#if 0
CPPUNIT_TEST( test3DSurfInterpP0P0PL_1 );
CPPUNIT_TEST( test3DSurfInterpP0P1_1 );
CPPUNIT_TEST( test3DSurfInterpP0P1PL_1 );
CPPUNIT_TEST( testInterpolationCU1D );
CPPUNIT_TEST( testInterpolationCU2D );
CPPUNIT_TEST( testInterpolationCU3D );
-#endif
CPPUNIT_TEST( test3DInterpP0P0_1 );
-#if 0 //dp
CPPUNIT_TEST( test3DInterpP0P0PL_1 );
CPPUNIT_TEST( test3DInterpP0P0PL_2 );
CPPUNIT_TEST( test3DInterpP0P0PL_3 );
CPPUNIT_TEST( test3DSurfInterpP1P0Bary_1 );
CPPUNIT_TEST( test3DInterpP1P0Bary_1 );
CPPUNIT_TEST( test3DTo1DInterpP0P0PL_1 );
-#endif
- CPPUNIT_TEST( test3D2DInterpP0P0_1 );
+ CPPUNIT_TEST( test3D2DBasicInterpP0P0 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_1 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_2 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_3 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_4 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_5 );
+ CPPUNIT_TEST( test3D2QuadHexaInterpP0P0_6 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_1 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_2 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_3 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_4 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_5 );
+ CPPUNIT_TEST( test3D2TriHexaInterpP0P0_6 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_1 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_2 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_3 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_4 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_5 );
+ CPPUNIT_TEST( test3D2QuadTetraInterpP0P0_6 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_1 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_2 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_3 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_4 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_5 );
+ CPPUNIT_TEST( test3D2TriTetraInterpP0P0_6 );
-#if 0//dp
CPPUNIT_TEST( test1DInterp_1 );
CPPUNIT_TEST( test2DCurveInterpP0P0_1 );
CPPUNIT_TEST( test2DCurveInterpP0P0_2 );
CPPUNIT_TEST( test2DCurveInterpP0P1_1 );
CPPUNIT_TEST( test2DCurveInterpP1P0_1 );
CPPUNIT_TEST( test2DCurveInterpP1P1_1 );
-#endif
CPPUNIT_TEST_SUITE_END();
public:
//MEDCouplingBasicsTest1.cxx
void test2DCurveInterpP1P0_1();
void test2DCurveInterpP1P1_1();
-#if 1//dp
- void test3D2DInterpP0P0_1();
-#endif
+ void test3D2DBasicInterpP0P0();
+ void test3D2QuadHexaInterpP0P0_1();
+ void test3D2QuadHexaInterpP0P0_2();
+ void test3D2QuadHexaInterpP0P0_3();
+ void test3D2QuadHexaInterpP0P0_4();
+ void test3D2QuadHexaInterpP0P0_5();
+ void test3D2QuadHexaInterpP0P0_6();
+ void test3D2TriHexaInterpP0P0_1();
+ void test3D2TriHexaInterpP0P0_2();
+ void test3D2TriHexaInterpP0P0_3();
+ void test3D2TriHexaInterpP0P0_4();
+ void test3D2TriHexaInterpP0P0_5();
+ void test3D2TriHexaInterpP0P0_6();
+ void test3D2QuadTetraInterpP0P0_1();
+ void test3D2QuadTetraInterpP0P0_2();
+ void test3D2QuadTetraInterpP0P0_3();
+ void test3D2QuadTetraInterpP0P0_4();
+ void test3D2QuadTetraInterpP0P0_5();
+ void test3D2QuadTetraInterpP0P0_6();
+ void test3D2TriTetraInterpP0P0_1();
+ void test3D2TriTetraInterpP0P0_2();
+ void test3D2TriTetraInterpP0P0_3();
+ void test3D2TriTetraInterpP0P0_4();
+ void test3D2TriTetraInterpP0P0_5();
+ void test3D2TriTetraInterpP0P0_6();
public:
static MEDCouplingUMesh *build3DSourceMesh_2();
static MEDCouplingUMesh *buildHexa8Mesh_1();
static MEDCouplingUMesh *buildPointe_1(MEDCouplingUMesh *&m1);
-#if 1//dp
- static MEDCouplingUMesh *build3D2DSourceMesh_1();
- static MEDCouplingUMesh *build3D2DTargetMesh_1();
-#endif
+ static MEDCouplingUMesh *build3D2DSourceMesh();
+ static MEDCouplingUMesh *build3D2DTargetMesh();
+ static MEDCouplingUMesh* build3D2DQuadSourceMesh(const double shiftX = 0.,
+ const double inclinationX = 0.);
+ static MEDCouplingUMesh* build3D2DTriSourceMesh(const double shiftX = 0.,
+ const double inclinationX = 0.);
+ static MEDCouplingUMesh* build3D2DTetraTargetMesh(const double inclinaisonX = 0.);
+ static MEDCouplingUMesh* build3D2DHexaTargetMesh(const double inclinaisonX = 0.);
static DataArrayDouble *buildCoordsForMultiTypes_1();
static MEDCouplingMultiFields *buildMultiFields_1();
static std::vector<MEDCouplingFieldDouble *> buildMultiFields_2();
static double sumAll(const std::vector< std::map<int,double> >& matrix);
+
+ private:
+ static int countNonZero(const std::vector< std::map<int,double> >& matrix);
+
+ static void test3D2DMeshesIntersection(MEDCouplingUMesh *sourceMesh,
+ MEDCouplingUMesh *targetMesh,
+ const double correctSurf,
+ const int correctDuplicateFacesNbr,
+ const int correctTotalIntersectFacesNbr = -1);
};
}
return ret;
}
-#if 1//dp
-
-#if 1
-MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1()
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh()
{
double sourceCoords[63]={-12., 6., 10., -12.,10., 6., -16.,10. , 10.,
-20., 0., 0., -12., 0., 0., -12., 0. , -4., -20.,0.,-4.,
return sourceMesh;
}
-MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh_1()
+MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh()
+{
+ double targetCoords[45]={-20., 0., 0., -20.,10., 0., -12.,10., 0.,
+ -12., 0., 0., -20., 0.,10., -20.,10.,10.,
+ -12.,10.,10., -12., 0.,10., -20., 0.,18.,
+ -20.,-5.,10., -20.,-5.,-4., -12.,-5.,-4.,
+ -12.,-5.,10., -20., 0.,-4., -12., 0.,-4.
+ };
+ int targetConn[20]={4,5,7,8, 0,3,2,1,4,7,6,5, 4,13,14,7,9,10,11,12};
+ MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+ targetMesh->setMeshDimension(3);
+ targetMesh->allocateCells(3);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,targetConn + 4);
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,targetConn + 12);
+ targetMesh->finishInsertingCells();
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(15,3);
+ std::copy(targetCoords,targetCoords+45,myCoords->getPointer());
+ targetMesh->setCoords(myCoords);
+ myCoords->decrRef();
+ return targetMesh;
+}
+
+MEDCouplingUMesh* MEDCouplingBasicsTest::build3D2DQuadSourceMesh(const double shiftX,
+ const double inclinationX)
{
- double sourceCoords[27]={-20.,0.,0., -20.,10.,0., -12.,10.,0., -12.,0.,0., -20.,0.,10., -20.,10.,10., -12.,10.,10., -12.,0.,10., -20.,0.,18.};
- int sourceConn[12]={4,5,7,8, 0,3,2,1,4,7,6,5};
MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New();
- sourceMesh->setMeshDimension(3);
- sourceMesh->allocateCells(2);
- sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,sourceConn);
- sourceMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,sourceConn + 4);
+ sourceMesh->setMeshDimension(2);
+
+ const int nbY = 4;
+ const int nbZ = 5;
+ const int nbYP1 = nbY + 1;
+ const int nbZP1 = nbZ + 1;
+ sourceMesh->allocateCells(nbY * nbZ);
+
+ int sourceConn[4];
+ for (int iY = 0; iY < nbY; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZ; ++iZ)
+ {
+ sourceConn[0] = iZ + iY * nbZP1;
+ sourceConn[1] = iZ + 1 + iY * nbZP1;
+ sourceConn[2] = iZ + 1 + (iY + 1) * nbZP1;
+ sourceConn[3] = iZ + (iY + 1) * nbZP1;
+ sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn);
+ }
+ }
sourceMesh->finishInsertingCells();
+
+ std::vector<double> sourceCoords;
+ for (int iY = 0; iY < nbYP1; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZP1; ++iZ)
+ {
+ sourceCoords.push_back(iY * inclinationX + shiftX);
+ sourceCoords.push_back(iY * 4.);
+ sourceCoords.push_back(iZ * 3.);
+ }
+
+ }
DataArrayDouble *myCoords=DataArrayDouble::New();
- myCoords->alloc(9,3);
- std::copy(sourceCoords,sourceCoords+27,myCoords->getPointer());
+ myCoords->alloc(nbYP1 * nbZP1,3);
+ std::copy(sourceCoords.begin(),sourceCoords.end(),myCoords->getPointer());
sourceMesh->setCoords(myCoords);
myCoords->decrRef();
+
return sourceMesh;
}
-#else // détruire dessous
+MEDCouplingUMesh* MEDCouplingBasicsTest::build3D2DTriSourceMesh(const double shiftX,
+ const double inclinationX)
+{
+ MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New();
+ sourceMesh->setMeshDimension(2);
+
+ const int nbY = 4;
+ const int nbZ = 5;
+ const int nbYP1 = nbY + 1;
+ const int nbZP1 = nbZ + 1;
+ sourceMesh->allocateCells(nbY * nbZ * 2);
+
+ int sourceConn[3];
+ for (int iY = 0; iY < nbY; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZ; ++iZ)
+ {
+ sourceConn[0] = iZ + iY * nbZP1;
+ sourceConn[1] = iZ + 1 + iY * nbZP1;
+ sourceConn[2] = iZ + 1 + (iY + 1) * nbZP1;
+ sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn);
+ sourceConn[0] = iZ + iY * nbZP1;
+ sourceConn[1] = iZ + (iY + 1) * nbZP1;
+ sourceConn[2] = iZ + 1 + (iY + 1) * nbZP1;
+ sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,sourceConn);
+ }
+ }
+ sourceMesh->finishInsertingCells();
+
+ std::vector<double> sourceCoords;
+ for (int iY = 0; iY < nbYP1; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZP1; ++iZ)
+ {
+ sourceCoords.push_back(iY * inclinationX + shiftX);
+ sourceCoords.push_back(iY * 4.);
+ sourceCoords.push_back(iZ * 3.);
+ }
-MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh_1()
+ }
+ DataArrayDouble *myCoords=DataArrayDouble::New();
+ myCoords->alloc(nbYP1 * nbZP1,3);
+ std::copy(sourceCoords.begin(),sourceCoords.end(),myCoords->getPointer());
+ sourceMesh->setCoords(myCoords);
+ myCoords->decrRef();
+
+ return sourceMesh;
+}
+
+MEDCouplingUMesh* MEDCouplingBasicsTest::build3D2DHexaTargetMesh(const double inclinationX)
{
- double targetCoords[12]={0.,0.,0, 1.,0.,0., 0.,1.,0., 0.,0.,1.};
- int targetConn[4]={0,1,2,3};
MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
targetMesh->setMeshDimension(3);
- targetMesh->allocateCells(1);
- targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+
+ const int nbX = 5;
+ const int nbY = 4;
+ const int nbZ = 5;
+ const int nbXP1 = nbX + 1;
+ const int nbYP1 = nbY + 1;
+ const int nbZP1 = nbZ + 1;
+ targetMesh->allocateCells(nbX * nbY * nbZ);
+
+ int targetConn[8];
+ for (int iX = 0; iX < nbX; ++iX)
+ {
+ for (int iY = 0; iY < nbY; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZ; ++iZ)
+ {
+ targetConn[0] = iZ + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[1] = iZ + 1 + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[2] = iZ + 1 + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetConn[3] = iZ + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetConn[4] = iZ + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[5] = iZ + 1 + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[6] = iZ + 1 + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[7] = iZ + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,targetConn);
+ }
+ }
+ }
targetMesh->finishInsertingCells();
+
+ std::vector<double> targetCoords;
+ for (int iX = 0; iX < nbXP1; ++iX)
+ {
+ for (int iY = 0; iY < nbYP1; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZP1; ++iZ)
+ {
+ targetCoords.push_back(iX * 3. + iY * inclinationX);
+ targetCoords.push_back(iY * 4.);
+ targetCoords.push_back(iZ * 3.);
+ }
+ }
+ }
DataArrayDouble *myCoords=DataArrayDouble::New();
- myCoords->alloc(4,3);
- std::copy(targetCoords,targetCoords+12,myCoords->getPointer());
+ myCoords->alloc(nbXP1 * nbYP1 * nbZP1, 3);
+ std::copy(targetCoords.begin(),targetCoords.end(),myCoords->getPointer());
targetMesh->setCoords(myCoords);
myCoords->decrRef();
+
return targetMesh;
}
-MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1()
+MEDCouplingUMesh* MEDCouplingBasicsTest::build3D2DTetraTargetMesh(const double inclinationX)
{
- double targetCoords[9]={0.5,0.,0., 0.5,0.5,0., 0.5,0.,0.5};
- int targetConn[3]={0,1,2};
MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
- targetMesh->setMeshDimension(2);
- targetMesh->allocateCells(1);
- targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn);
+ targetMesh->setMeshDimension(3);
+
+ const int nbX = 5;
+ const int nbY = 4;
+ const int nbZ = 5;
+ const int nbXP1 = nbX + 1;
+ const int nbYP1 = nbY + 1;
+ const int nbZP1 = nbZ + 1;
+ targetMesh->allocateCells(nbX * nbY * nbZ * 5);
+
+ int targetConn[4];
+ for (int iX = 0; iX < nbX; ++iX)
+ {
+ for (int iY = 0; iY < nbY; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZ; ++iZ)
+ {
+ targetConn[0] = iZ + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[1] = iZ + 1 + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[2] = iZ + 1 + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[3] = iZ + 1 + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ targetConn[0] = iZ + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[1] = iZ + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[2] = iZ + 1 + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[3] = iZ + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ targetConn[0] = iZ + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[1] = iZ + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetConn[2] = iZ + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[3] = iZ + 1 + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ targetConn[0] = iZ + 1 + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[1] = iZ + 1 + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[2] = iZ + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[3] = iZ + 1 + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ targetConn[0] = iZ + ( iY + iX * nbYP1) * nbZP1;
+ targetConn[1] = iZ + 1 + ((iY + 1) + iX * nbYP1) * nbZP1;
+ targetConn[2] = iZ + 1 + ( iY + (iX + 1) * nbYP1) * nbZP1;
+ targetConn[3] = iZ + ((iY + 1) + (iX + 1) * nbYP1) * nbZP1;
+ targetMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,targetConn);
+ }
+ }
+ }
targetMesh->finishInsertingCells();
+
+ std::vector<double> targetCoords;
+ for (int iX = 0; iX < nbXP1; ++iX)
+ {
+ for (int iY = 0; iY < nbYP1; ++iY)
+ {
+ for (int iZ = 0; iZ < nbZP1; ++iZ)
+ {
+ targetCoords.push_back(iX * 3. + iY * inclinationX);
+ targetCoords.push_back(iY * 4.);
+ targetCoords.push_back(iZ * 3.);
+ }
+ }
+ }
DataArrayDouble *myCoords=DataArrayDouble::New();
- myCoords->alloc(3,3);
- std::copy(targetCoords,targetCoords+9,myCoords->getPointer());
+ myCoords->alloc(nbXP1 * nbYP1 * nbZP1, 3);
+ std::copy(targetCoords.begin(),targetCoords.end(),myCoords->getPointer());
targetMesh->setCoords(myCoords);
myCoords->decrRef();
+
return targetMesh;
}
-#endif
-#endif
using namespace ParaMEDMEM;
+typedef std::vector<std::map<int,double> > IntersectionMatrix;
+
void MEDCouplingBasicsTest::test2DInterpP0P0_1()
{
MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1();
targetMesh->decrRef();
}
-#if 1//dp
-#if 0
-void MEDCouplingBasicsTest::test3D2DInterpP0P0_1()
+void MEDCouplingBasicsTest::test3D2DBasicInterpP0P0()
{
- MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh_1();
- MEDCouplingUMesh *targetMesh=build3D2DTargetMesh_1();
+ MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh();
+ MEDCouplingUMesh *targetMesh=build3D2DTargetMesh();
MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
INTERP_KERNEL::Interpolation3D2D myInterpolator;
myInterpolator.setPrecision(1e-12);
- std::vector<std::map<int,double> > res;
+ std::vector<std::map<int,double> > matrix;
INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 };
- for ( int i = 0; i < 4; ++i )
+ for ( size_t i = 0; i < sizeof(sp)/sizeof(sp[0]); ++i )
{
myInterpolator.setSplittingPolicy( sp[i] );
- res.clear();
- myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
- CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][0],1e-12);
- //CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7);
+ matrix.clear();
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,matrix,"P0P0");
+
+ CPPUNIT_ASSERT_EQUAL(3,(int)matrix.size());
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[0][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[0][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,matrix[0][2],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(8. ,matrix[0][3],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(2.5 ,matrix[0][4],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[0][5],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(32. ,matrix[0][6],1e-12);
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(8.*sqrt(3.),matrix[1][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[1][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,matrix[1][2],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,matrix[1][3],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[1][4],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,matrix[1][5],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,matrix[1][6],1e-12);
+
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[2][0],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(32. ,matrix[2][1],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[2][2],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[2][3],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,matrix[2][4],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,matrix[2][5],1e-12);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(112. ,matrix[2][6],1e-12);
+
+ INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces = myInterpolator.retrieveDuplicateFaces();
+ CPPUNIT_ASSERT_EQUAL(3,(int)duplicateFaces.size());
+
+ INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType correctDuplicateFaces;
+ std::set<int> face2;
+ face2.insert(0);
+ face2.insert(1);
+ correctDuplicateFaces[2] = face2;
+ std::set<int> face5;
+ face5.insert(1);
+ face5.insert(2);
+ correctDuplicateFaces[5] = face5;
+ std::set<int> face6;
+ face6.insert(0);
+ face6.insert(1);
+ face6.insert(2);
+ correctDuplicateFaces[6] = face6;
+
+ CPPUNIT_ASSERT(correctDuplicateFaces == duplicateFaces);
}
//clean up
sourceMesh->decrRef();
targetMesh->decrRef();
}
-#else
-void MEDCouplingBasicsTest::test3D2DInterpP0P0_1()
+
+int MEDCouplingBasicsTest::countNonZero(const std::vector< std::map<int,double> >& matrix)
{
- MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh_1();
- MEDCouplingUMesh *targetMesh=build3D2DTargetMesh_1();
+ int ret=0.;
+ for(std::vector< std::map<int,double> >::const_iterator iter=matrix.begin();iter!=matrix.end();iter++)
+ for(std::map<int,double>::const_iterator iter2=(*iter).begin();iter2!=(*iter).end();iter2++)
+ if (!INTERP_KERNEL::epsilonEqual((*iter2).second, 0.)) ret +=1;
+ return ret;
+}
+void MEDCouplingBasicsTest::test3D2DMeshesIntersection(MEDCouplingUMesh *sourceMesh,
+ MEDCouplingUMesh *targetMesh,
+ const double correctSurf,
+ const int correctDuplicateFacesNbr,
+ const int correctTotalIntersectFacesNbr)
+{
MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(sourceMesh);
MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
INTERP_KERNEL::Interpolation3D2D myInterpolator;
myInterpolator.setPrecision(1e-12);
- std::vector<std::map<int,double> > res;
- //INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::GENERAL_48 };
+ const double prec = 1.0e-5;
+ IntersectionMatrix matrix;
INTERP_KERNEL::SplittingPolicy sp[] = { INTERP_KERNEL::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 };
- for ( int i = 0; i < 1; ++i )
+ for ( size_t i = 0; i < sizeof(sp)/sizeof(sp[0]); ++i )
{
myInterpolator.setSplittingPolicy( sp[i] );
- res.clear();
- myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0");
+ matrix.clear();
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,matrix,"P0P0");
- CPPUNIT_ASSERT_EQUAL(2,(int)res.size());
-
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][0],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][1],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[0][2],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(8. ,res[0][3],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(2.5 ,res[0][4],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[0][5],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(32. ,res[0][6],1e-12);
-
- CPPUNIT_ASSERT_DOUBLES_EQUAL(8.*sqrt(3.),res[1][0],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[1][1],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(40. ,res[1][2],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][3],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0. ,res[1][4],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][5],1e-12);
- CPPUNIT_ASSERT_DOUBLES_EQUAL(80. ,res[1][6],1e-12);
-
- //CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7);
+ std::cout.precision(16);
+
+ const double surf = sumAll(matrix);
+ LOG(1, "surf = " << surf <<" correctSurf = " << correctSurf );
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(correctSurf, surf, prec * std::max(correctSurf, surf));
+
+ INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces = myInterpolator.retrieveDuplicateFaces();
+ int duplicateFacesNbr = duplicateFaces.size();
+ LOG(1, "duplicateFacesNbr = " << duplicateFacesNbr <<" correctDuplicateFacesNbr = " << correctDuplicateFacesNbr);
+ CPPUNIT_ASSERT_EQUAL(correctDuplicateFacesNbr, duplicateFacesNbr);
+
+ if (correctTotalIntersectFacesNbr >= 0)
+ {
+ int totalIntersectFacesNbr = countNonZero(matrix);
+ LOG(1, "totalIntersectFacesNbr = " << totalIntersectFacesNbr <<" correctTotalIntersectFacesNbr = " << correctTotalIntersectFacesNbr );
+ CPPUNIT_ASSERT_EQUAL(correctTotalIntersectFacesNbr, totalIntersectFacesNbr);
+ }
}
//clean up
sourceMesh->decrRef();
targetMesh->decrRef();
}
-#endif
-#endif
+
+//dp TODO DP : adapter les commentaires
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh();
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 240., 0, 20);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_2()
+{
+ const double shiftX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 240., 20, 2 * 20);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_3()
+{
+ const double shiftX = 1.5;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 300., 0, 20);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_4()
+{
+ const double shiftX = 3.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 300., 20, 2 * 20);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_5()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 180., 0, 15);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadHexaInterpP0P0_6()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 150., 0, 10);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh();
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 240., 0, 40);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_2()
+{
+ const double shiftX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 240., 40, 2 * 40);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_3()
+{
+ const double shiftX = 1.5;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 300., 0, 40);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_4()
+{
+ const double shiftX = 3.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 300., 40, 2 * 40);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_5()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 180., 0, 30);
+}
+
+void MEDCouplingBasicsTest::test3D2TriHexaInterpP0P0_6()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DHexaTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 150., 0, 20);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh();
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 240., 20, 40);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_2()
+{
+ const double shiftX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 240., 20, 2 * 40);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_3()
+{
+ const double shiftX = 1.5;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 300., 0, 100);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_4()
+{
+ const double shiftX = 3.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 300., 20, 2 * 40);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_5()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 180., 0, 45);
+}
+
+void MEDCouplingBasicsTest::test3D2QuadTetraInterpP0P0_6()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DQuadSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 150., 0, 30);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_1()
+{
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh();
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 240., 0, 40);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_2()
+{
+ const double shiftX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 240., 40, 40 + 80);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_3()
+{
+ const double shiftX = 1.5;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 300., 0);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_4()
+{
+ const double shiftX = 3.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 2. * 300., 40, 40 + 80);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_5()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh(inclinationX);
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 180., 0);
+}
+
+void MEDCouplingBasicsTest::test3D2TriTetraInterpP0P0_6()
+{
+ const double shiftX = 9.;
+ const double inclinationX = 3.;
+ MEDCouplingUMesh *sourceMesh=build3D2DTriSourceMesh(shiftX, inclinationX);
+ MEDCouplingUMesh *targetMesh=build3D2DTetraTargetMesh();
+ test3D2DMeshesIntersection(sourceMesh, targetMesh, 150., 0);
+}
