From 48d2ff65df3dc1a888f8ce377ad254a1b6b4ff73 Mon Sep 17 00:00:00 2001 From: bph Date: Fri, 5 Aug 2011 14:46:18 +0000 Subject: [PATCH] *** empty log message *** --- src/INTERP_KERNEL/Interpolation.hxx | 2 + src/INTERP_KERNEL/Interpolation.txx | 39 ++ src/INTERP_KERNEL/Interpolation3D2D.hxx | 17 +- src/INTERP_KERNEL/Interpolation3D2D.txx | 52 +-- src/INTERP_KERNEL/InterpolationPlanar.hxx | 1 + src/INTERP_KERNEL/InterpolationPlanar.txx | 26 +- src/INTERP_KERNEL/InterpolationUtils.hxx | 3 - .../Polyhedron3D2DIntersectorP0P0.hxx | 24 +- .../Polyhedron3D2DIntersectorP0P0.txx | 89 ++++- src/INTERP_KERNEL/SplitterTetra.hxx | 24 +- src/INTERP_KERNEL/SplitterTetra.txx | 187 ++++------ src/INTERP_KERNEL/TransformedTriangle.cxx | 9 +- src/INTERP_KERNEL/TransformedTriangle.hxx | 2 - src/INTERP_KERNEL/VectorUtils.hxx | 2 - src/INTERP_KERNELTest/MeshTestToolkit.hxx | 2 +- src/INTERP_KERNELTest/TestInterpKernel.cxx | 2 + .../Test/MEDCouplingBasicsTest.hxx | 83 ++++- .../Test/MEDCouplingBasicsTest0.cxx | 250 +++++++++++-- .../Test/MEDCouplingBasicsTestInterp.cxx | 339 +++++++++++++++--- 19 files changed, 854 insertions(+), 299 deletions(-) diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx index 128ca6e3c..4a34b863a 100644 --- a/src/INTERP_KERNEL/Interpolation.hxx +++ b/src/INTERP_KERNEL/Interpolation.hxx @@ -43,6 +43,8 @@ namespace INTERP_KERNEL template 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 + static double CalculateCharacteristicSizeOfMeshes(const MyMeshType& myMeshS, const MyMeshType& myMeshT, const int printLevel); protected: template int fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method); diff --git a/src/INTERP_KERNEL/Interpolation.txx b/src/INTERP_KERNEL/Interpolation.txx index fc5d4c30e..a5bdca7e4 100644 --- a/src/INTERP_KERNEL/Interpolation.txx +++ b/src/INTERP_KERNEL/Interpolation.txx @@ -22,6 +22,7 @@ #include "Interpolation.hxx" #include "IntegralUniformIntersector.hxx" #include "IntegralUniformIntersector.txx" +#include "VectorUtils.hxx" namespace INTERP_KERNEL { @@ -70,6 +71,44 @@ namespace INTERP_KERNEL srcMeth=methodC.substr(0,2); trgMeth=methodC.substr(2); } + + template + template + double Interpolation::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::max(); + if(nbMailleS!=0) + { + diagonalS=getDistanceBtw2Pts(BoxS+SPACEDIM,BoxS); + dimCaracteristicS=diagonalS/nbMailleS; + } + double diagonalT,dimCaracteristicT=std::numeric_limits::max(); + if(nbMailleT!=0) + { + diagonalT=getDistanceBtw2Pts(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 diff --git a/src/INTERP_KERNEL/Interpolation3D2D.hxx b/src/INTERP_KERNEL/Interpolation3D2D.hxx index 3eabde9b5..f38aab642 100644 --- a/src/INTERP_KERNEL/Interpolation3D2D.hxx +++ b/src/INTERP_KERNEL/Interpolation3D2D.hxx @@ -20,6 +20,9 @@ #ifndef __INTERPOLATION3D2D_HXX__ #define __INTERPOLATION3D2D_HXX__ +#include +#include + #include "Interpolation.hxx" #include "NormalizedUnstructuredMesh.hxx" #include "InterpolationOptions.hxx" @@ -29,12 +32,22 @@ namespace INTERP_KERNEL class Interpolation3D2D : public Interpolation { public: + typedef typename std::map > DuplicateFacesType; + Interpolation3D2D(); Interpolation3D2D(const InterpolationOptions& io); - template - int interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method); + template + 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; }; } diff --git a/src/INTERP_KERNEL/Interpolation3D2D.txx b/src/INTERP_KERNEL/Interpolation3D2D.txx index 8bfd9f071..27ddb7e8f 100644 --- a/src/INTERP_KERNEL/Interpolation3D2D.txx +++ b/src/INTERP_KERNEL/Interpolation3D2D.txx @@ -33,20 +33,9 @@ #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 - -#else // use BBTree class #include "BBTree.txx" -#endif - namespace INTERP_KERNEL { /** @@ -67,11 +56,14 @@ 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 - int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method) + template + 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 @@ -84,6 +76,7 @@ namespace INTERP_KERNEL std::vector*> targetElems(numTargetElems); std::map*, int> indices; + DuplicateFacesType intersectFaces; for(unsigned long i = 0 ; i < numSrcElems ; ++i) srcElems[i] = new MeshElement(i, srcMesh); @@ -91,15 +84,20 @@ namespace INTERP_KERNEL for(unsigned long i = 0 ; i < numTargetElems ; ++i) targetElems[i] = new MeshElement(i, targetMesh); - Intersector3D* intersector=0; + Intersector3D* 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(targetMesh, srcMesh, - getSplittingPolicy()); + intersector=new Polyhedron3D2DIntersectorP0P0(targetMesh, + srcMesh, + dimCaracteristic, + getPrecision(), + intersectFaces, + getSplittingPolicy()); break; default: throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P0P0 interp specified : must be Triangulation."); @@ -108,9 +106,7 @@ namespace INTERP_KERNEL 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 @@ -131,11 +127,7 @@ namespace INTERP_KERNEL 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 @@ -158,12 +150,22 @@ namespace INTERP_KERNEL 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(); diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx index 5699fbc48..c47a0205e 100755 --- a/src/INTERP_KERNEL/InterpolationPlanar.hxx +++ b/src/INTERP_KERNEL/InterpolationPlanar.hxx @@ -43,6 +43,7 @@ namespace INTERP_KERNEL // Main function to interpolate triangular and quadratic meshes template 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(); } diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx index a8ea915f7..23df2913e 100644 --- a/src/INTERP_KERNEL/InterpolationPlanar.txx +++ b/src/INTERP_KERNEL/InterpolationPlanar.txx @@ -86,8 +86,7 @@ namespace INTERP_KERNEL 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 @@ -123,32 +122,15 @@ namespace INTERP_KERNEL /***********************************************************/ 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::max(); - if(nbMailleS!=0) - { - diagonalS=getDistanceBtw2Pts(BoxS+SPACEDIM,BoxS); - dimCaracteristicS=diagonalS/nbMailleS; - } - double diagonalT,dimCaracteristicT=std::numeric_limits::max(); - if(nbMailleT!=0) - { - diagonalT=getDistanceBtw2Pts(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; } diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx index 33c035227..c79b0ce9e 100644 --- a/src/INTERP_KERNEL/InterpolationUtils.hxx +++ b/src/INTERP_KERNEL/InterpolationUtils.hxx @@ -22,9 +22,6 @@ #include "INTERPKERNELDefines.hxx" #include "InterpKernelException.hxx" -#if 1//dp -#include "VectorUtils.hxx" -#endif #include "NormalizedUnstructuredMesh.hxx" diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx index 526493e6e..c94c51cb0 100644 --- a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx @@ -33,22 +33,31 @@ namespace INTERP_KERNEL * the source elements. * */ - template - class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0 - { + template + class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0 + { + typedef typename std::map > 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& srcCells, MyMatrix& res); + void intersectCells(ConnType targetCell, + const std::vector& srcCells, + MyMatrixType& matrix); private: void releaseArrays(); @@ -59,6 +68,11 @@ namespace INTERP_KERNEL SplitterTetra2 _split; + double _dim_caracteristic; + double _precision; + + DuplicateFacesType& _intersect_faces; + }; } diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx index ecddf3c6e..12a17bd45 100644 --- a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx @@ -38,10 +38,18 @@ namespace INTERP_KERNEL * @param srcMesh mesh containing the source elements * @param policy splitting policy to be used */ - template - Polyhedron3D2DIntersectorP0P0::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, const MyMeshType& srcMesh, - SplittingPolicy policy) - : Intersector3DP0P0(targetMesh,srcMesh),_split(targetMesh,srcMesh,policy) + template + Polyhedron3D2DIntersectorP0P0::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, + const MyMeshType& srcMesh, + const double dimCaracteristic, + const double precision, + DuplicateFacesType& intersectFaces, + SplittingPolicy policy) + : Intersector3DP0P0(targetMesh,srcMesh), + _split(targetMesh,srcMesh,policy), + _dim_caracteristic(dimCaracteristic), + _precision(precision), + _intersect_faces(intersectFaces) { } @@ -50,14 +58,14 @@ namespace INTERP_KERNEL * Liberates the SplitterTetra objects and potential sub-node points that have been allocated. * */ - template - Polyhedron3D2DIntersectorP0P0::~Polyhedron3D2DIntersectorP0P0() + template + Polyhedron3D2DIntersectorP0P0::~Polyhedron3D2DIntersectorP0P0() { releaseArrays(); } - template - void Polyhedron3D2DIntersectorP0P0::releaseArrays() + template + void Polyhedron3D2DIntersectorP0P0::releaseArrays() { for(typename std::vector< SplitterTetra* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) delete *iter; @@ -75,25 +83,28 @@ namespace INTERP_KERNEL * @param srcCells in C mode. * */ - template - void Polyhedron3D2DIntersectorP0P0::intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res) + template + void Polyhedron3D2DIntersectorP0P0::intersectCells(ConnType targetCell, + const std::vector& srcCells, + MyMatrixType& matrix) { - int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)); + int nbOfNodesT=Intersector3D::_target_mesh.getNumberOfNodesOfElement(OTT::indFC(targetCell)); releaseArrays(); _split.splitTargetCell(targetCell,nbOfNodesT,_tetra); for(typename std::vector::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++) { double surface = 0.; - std::set listOfTetraFacesTreated; + std::multiset listOfTetraFacesTreated; + std::set listOfTetraFacesColinear; // calculate the coordinates of the nodes const NumberingPolicy numPol=MyMeshType::My_numPol; typename MyMeshType::MyConnType cellSrc = *iterCellS; int cellSrcIdx = OTT::indFC(cellSrc); - NormalizedCellType normCellType=Intersector3D::_src_mesh.getTypeOfElement(cellSrcIdx); + NormalizedCellType normCellType=Intersector3D::_src_mesh.getTypeOfElement(cellSrcIdx); const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); - const MyMeshType& _src_mesh = Intersector3D::_src_mesh; + const MyMeshType& _src_mesh = Intersector3D::_src_mesh; unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(cellSrcIdx) : cellModelCell.getNumberOfNodes(); int *polyNodes=new int[nbOfNodes4Type]; double **polyCoords = new double*[nbOfNodes4Type]; @@ -106,10 +117,52 @@ namespace INTERP_KERNEL } for(typename std::vector*>::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::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::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 targetCellSet; + targetCellSet.insert(targetCell); + _intersect_faces.insert(std::make_pair(cellSrcIdx, targetCellSet)); + } + + } + + } + delete[] polyNodes; delete[] polyCoords; diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index e1bb11d0d..6e2191ff7 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -29,6 +29,8 @@ #include #include #include +#include +//dp#include namespace INTERP_KERNEL { @@ -194,7 +196,10 @@ namespace INTERP_KERNEL const int polyNodesNbr, const int *const polyNodes, const double *const *const polyCoords, - std::set& listOfTetraFacesTreated); + const double dimCaracteristic, + const double precision, + std::multiset& listOfTetraFacesTreated, + std::set& listOfTetraFacesColinear); double intersectTetra(const double** tetraCorners); @@ -210,14 +215,16 @@ namespace INTERP_KERNEL // 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); @@ -286,8 +293,6 @@ namespace INTERP_KERNEL isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); } -#if 1//dp - //dp à supprimer template inline void SplitterTetra::checkIsStrictlyOutside(const double* pt, bool* isStrictlyOutside, const double errTol) const { @@ -300,7 +305,6 @@ namespace INTERP_KERNEL 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. @@ -321,6 +325,16 @@ namespace INTERP_KERNEL _nodes[globalNodeNum] = transformedNode; } + //dp TODO DP : supprimer la précédente et commentariser + template + inline void SplitterTetra::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. diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index c3748f506..9e6659907 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -259,7 +259,6 @@ namespace INTERP_KERNEL faceNodes=new int[faceModel.getNumberOfNodes()]; cellModelCell.fillSonCellNodalConnectivity(ii,cellNodes,faceNodes); } - // intersect a son with the unit tetra switch(faceType) { @@ -375,9 +374,10 @@ namespace INTERP_KERNEL return std::fabs(1.0 / _t->determinant() * totalVolume) ; } -#if 1//dp TODO DP : à commenter + //dp TODO DP : à commenter template double SplitterTetra::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) @@ -479,57 +479,51 @@ namespace INTERP_KERNEL 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(&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 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 - bool SplitterTetra::isFacesCoplanar(const double *const plane_normal, const double plane_constant, + bool SplitterTetra::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. @@ -549,9 +543,6 @@ namespace INTERP_KERNEL 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. @@ -571,65 +562,15 @@ namespace INTERP_KERNEL const int polyNodesNbr, const int *const polyNodes, const double *const *const polyCoords, - std::set& listOfTetraFacesTreated) + const double dimCaracteristic, + const double precision, + std::multiset& listOfTetraFacesTreated, + std::set& 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 inter; - INTERP_KERNEL::intersec_de_triangle(coordsPoly[0], coordsPoly[1 + iTri], coordsPoly[2 + iTri], - coordsTriNode1, coordsTriNode2, coordsTriNode3, - inter, - 1., //dp inter, PlanarIntersector::_dim_caracteristic, - DEFAULT_ABS_TOL); //dp PlanarIntersector::_precision); - ConnType nb_inter=((ConnType)inter.size())/2; - if(nb_inter >3) inter=reconstruct_polygon(inter); - for(ConnType i = 1; i(&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) { @@ -645,24 +586,15 @@ namespace INTERP_KERNEL 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::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); @@ -704,7 +636,6 @@ namespace INTERP_KERNEL _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; @@ -731,18 +662,25 @@ namespace INTERP_KERNEL 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::_dim_caracteristic, - DEFAULT_ABS_TOL); //dp PlanarIntersector::_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 @@ -760,17 +698,17 @@ namespace INTERP_KERNEL { 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 : // @@ -793,12 +731,12 @@ namespace INTERP_KERNEL { 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 @@ -807,12 +745,12 @@ namespace INTERP_KERNEL { 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; @@ -828,11 +766,11 @@ namespace INTERP_KERNEL 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]; } } } @@ -850,14 +788,14 @@ namespace INTERP_KERNEL // 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; } /** @@ -1084,7 +1022,12 @@ namespace INTERP_KERNEL 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* t = new SplitterTetra(_src_mesh, nodes,conn); tetra.push_back(t); @@ -1130,12 +1073,8 @@ namespace INTERP_KERNEL 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* t = new SplitterTetra(_src_mesh, nodes,conn); tetra.push_back(t); @@ -1193,18 +1132,22 @@ namespace INTERP_KERNEL 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* t = new SplitterTetra(_src_mesh, nodes, conn); tetra.push_back(t); } diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 538b13a18..c0aec59a3 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -19,9 +19,8 @@ #include "TransformedTriangle.hxx" #include "VectorUtils.hxx" -#if 1//dp #include "TetraAffineTransform.hxx" -#endif +//dp #include "InterpolationUtils.hxx" #include #include #include @@ -589,7 +588,7 @@ namespace INTERP_KERNEL 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); @@ -797,6 +796,10 @@ namespace INTERP_KERNEL }; 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; } diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 2dba44370..a396b2beb 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -46,9 +46,7 @@ namespace INTERP_TEST 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 diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index 1fcf93e0e..3866913ef 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -80,8 +80,6 @@ namespace INTERP_KERNEL return ss.str(); } - - //dp : à supprimer /** * Subtracts two a double[3] - vectors and store the result in res * diff --git a/src/INTERP_KERNELTest/MeshTestToolkit.hxx b/src/INTERP_KERNELTest/MeshTestToolkit.hxx index e6db9bcda..61d0db78a 100644 --- a/src/INTERP_KERNELTest/MeshTestToolkit.hxx +++ b/src/INTERP_KERNELTest/MeshTestToolkit.hxx @@ -39,7 +39,7 @@ namespace INTERP_KERNEL namespace MEDMEM { class MESH; -}; +} namespace INTERP_TEST { diff --git a/src/INTERP_KERNELTest/TestInterpKernel.cxx b/src/INTERP_KERNELTest/TestInterpKernel.cxx index 25689c8c4..45a2772c6 100644 --- a/src/INTERP_KERNELTest/TestInterpKernel.cxx +++ b/src/INTERP_KERNELTest/TestInterpKernel.cxx @@ -25,6 +25,7 @@ #include "TransformedTriangleIntersectTest.hxx" #include "TransformedTriangleTest.hxx" #include "UnitTetraIntersectionBaryTest.hxx" +#include "UnitTetra3D2DIntersectionTest.hxx" #ifdef DISABLE_MICROMED #include "HexaTests.hxx" @@ -46,6 +47,7 @@ CPPUNIT_TEST_SUITE_REGISTRATION( SingleElementPlanarTests ); 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 ); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 7b27d5f5f..769e88ccb 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -36,7 +36,6 @@ namespace ParaMEDMEM { CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest); //MEDCouplingBasicsTest1.cxx -#if 0//dp CPPUNIT_TEST( testArray ); CPPUNIT_TEST( testArray2 ); CPPUNIT_TEST( testArray3 ); @@ -238,9 +237,7 @@ namespace ParaMEDMEM 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 ); @@ -255,10 +252,8 @@ namespace ParaMEDMEM 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 ); @@ -279,18 +274,39 @@ namespace ParaMEDMEM 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 @@ -541,9 +557,31 @@ namespace ParaMEDMEM 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(); @@ -585,15 +623,28 @@ namespace ParaMEDMEM 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 buildMultiFields_2(); static double sumAll(const std::vector< std::map >& matrix); + + private: + static int countNonZero(const std::vector< std::map >& matrix); + + static void test3D2DMeshesIntersection(MEDCouplingUMesh *sourceMesh, + MEDCouplingUMesh *targetMesh, + const double correctSurf, + const int correctDuplicateFacesNbr, + const int correctTotalIntersectFacesNbr = -1); }; } diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx index ad76c84f3..c489a82da 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx @@ -1174,10 +1174,7 @@ double MEDCouplingBasicsTest::sumAll(const std::vector< std::map >& 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., @@ -1206,58 +1203,249 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1() 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 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 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 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 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 diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx index 16b053689..a80d9e9de 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -39,6 +39,8 @@ using namespace ParaMEDMEM; +typedef std::vector > IntersectionMatrix; + void MEDCouplingBasicsTest::test2DInterpP0P0_1() { MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1(); @@ -2267,74 +2269,327 @@ void MEDCouplingBasicsTest::test2DCurveInterpP1P1_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 > res; + std::vector > 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 face2; + face2.insert(0); + face2.insert(1); + correctDuplicateFaces[2] = face2; + std::set face5; + face5.insert(1); + face5.insert(2); + correctDuplicateFaces[5] = face5; + std::set 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 >& matrix) { - MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh_1(); - MEDCouplingUMesh *targetMesh=build3D2DTargetMesh_1(); + int ret=0.; + for(std::vector< std::map >::const_iterator iter=matrix.begin();iter!=matrix.end();iter++) + for(std::map::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 > 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); +} -- 2.39.2