From: bph Date: Tue, 12 Jul 2011 14:39:23 +0000 (+0000) Subject: lot 1.1 X-Git-Tag: EndWorkDidier~8 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=8281288a661b525ae02c8ea1ae944dc45f246b09;p=tools%2Fmedcoupling.git lot 1.1 --- diff --git a/src/INTERP_KERNEL/Interpolation3D2D.cxx b/src/INTERP_KERNEL/Interpolation3D2D.cxx new file mode 100644 index 000000000..d29392342 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.cxx @@ -0,0 +1,40 @@ +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#include "Interpolation3D2D.hxx" + +namespace INTERP_KERNEL +{ + /** + * \defgroup interpolation3D Interpolation3D + * \class Interpolation3D + * \brief Class used to calculate the volumes of intersection between the elements of two 3D meshes. + * + */ + /** + * Default constructor + * + */ + Interpolation3D2D::Interpolation3D2D() + { + } + Interpolation3D2D::Interpolation3D2D(const InterpolationOptions& io):Interpolation(io) + { + } +} diff --git a/src/INTERP_KERNEL/Interpolation3D2D.hxx b/src/INTERP_KERNEL/Interpolation3D2D.hxx new file mode 100644 index 000000000..3eabde9b5 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.hxx @@ -0,0 +1,41 @@ +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __INTERPOLATION3D2D_HXX__ +#define __INTERPOLATION3D2D_HXX__ + +#include "Interpolation.hxx" +#include "NormalizedUnstructuredMesh.hxx" +#include "InterpolationOptions.hxx" + +namespace INTERP_KERNEL +{ + class Interpolation3D2D : public Interpolation + { + public: + Interpolation3D2D(); + Interpolation3D2D(const InterpolationOptions& io); + template + int interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method); + private: + SplittingPolicy _splitting_policy; + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Interpolation3D2D.txx b/src/INTERP_KERNEL/Interpolation3D2D.txx new file mode 100644 index 000000000..8bfd9f071 --- /dev/null +++ b/src/INTERP_KERNEL/Interpolation3D2D.txx @@ -0,0 +1,185 @@ +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef __INTERPOLATION3D2D_TXX__ +#define __INTERPOLATION3D2D_TXX__ + +#include "Interpolation3D2D.hxx" +#include "Interpolation.txx" +#include "MeshElement.txx" +#include "TransformedTriangle.hxx" +#include "Polyhedron3D2DIntersectorP0P0.txx" +#include "PointLocator3DIntersectorP0P0.txx" +#include "PolyhedronIntersectorP0P1.txx" +#include "PointLocator3DIntersectorP0P1.txx" +#include "PolyhedronIntersectorP1P0.txx" +#include "PolyhedronIntersectorP1P0Bary.txx" +#include "PointLocator3DIntersectorP1P0.txx" +#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 +{ + /** + * Calculates the matrix of volumes of intersection between the elements of srcMesh and the elements of targetMesh. + * The calculation is done in two steps. First a filtering process reduces the number of pairs of elements for which the + * calculation must be carried out by eliminating pairs that do not intersect based on their bounding boxes. Then, the + * volume of intersection is calculated by an object of type Intersector3D for the remaining pairs, and entered into the + * intersection matrix. + * + * The matrix is partially sparse : it is a vector of maps of integer - double pairs. + * It can also be an INTERP_KERNEL::Matrix object. + * The length of the vector is equal to the number of target elements - for each target element there is a map, regardless + * of whether the element intersects any source elements or not. But in the maps there are only entries for those source elements + * which have a non-zero intersection volume with the target element. The vector has indices running from + * 0 to (nb target elements - 1), meaning that the map for target element i is stored at index i - 1. In the maps, however, + * the indexing is more natural : the intersection volume of the target element i with source element j is found at matrix[i-1][j]. + * + + * @param srcMesh 3-dimensional source mesh + * @param targetMesh 3-dimesional target mesh, containing only tetraedra + * @param result matrix in which the result is stored + * + */ + template + int Interpolation3D2D::interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method) + { + typedef typename MyMeshType::MyConnType ConnType; + // create MeshElement objects corresponding to each element of the two meshes + const unsigned long numSrcElems = srcMesh.getNumberOfElements(); + const unsigned long numTargetElems = targetMesh.getNumberOfElements(); + + LOG(2, "Source mesh has " << numSrcElems << " elements and target mesh has " << numTargetElems << " elements "); + + std::vector*> srcElems(numSrcElems); + std::vector*> targetElems(numTargetElems); + + std::map*, int> indices; + + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + srcElems[i] = new MeshElement(i, srcMesh); + + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + targetElems[i] = new MeshElement(i, targetMesh); + + Intersector3D* intersector=0; + std::string methC = InterpolationOptions::filterInterpolationMethod(method); + if(methC=="P0P0") + { + switch(InterpolationOptions::getIntersectionType()) + { + case Triangulation: + intersector=new Polyhedron3D2DIntersectorP0P0(targetMesh, srcMesh, + 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 + + // create BBTree structure + // - get bounding boxes + double* bboxes = new double[6 * numSrcElems]; + int* srcElemIdx = new int[numSrcElems]; + for(unsigned long i = 0; i < numSrcElems ; ++i) + { + // get source bboxes in right order + const BoundingBox* box = srcElems[i]->getBoundingBox(); + bboxes[6*i+0] = box->getCoordinate(BoundingBox::XMIN); + bboxes[6*i+1] = box->getCoordinate(BoundingBox::XMAX); + bboxes[6*i+2] = box->getCoordinate(BoundingBox::YMIN); + bboxes[6*i+3] = box->getCoordinate(BoundingBox::YMAX); + bboxes[6*i+4] = box->getCoordinate(BoundingBox::ZMIN); + bboxes[6*i+5] = box->getCoordinate(BoundingBox::ZMAX); + + // source indices have to begin with zero for BBox, I think + srcElemIdx[i] = srcElems[i]->getIndex(); + } + +#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 + for(unsigned long i = 0; i < numTargetElems; ++i) + { + const BoundingBox* box = targetElems[i]->getBoundingBox(); + const int targetIdx = targetElems[i]->getIndex(); + + // get target bbox in right order + double targetBox[6]; + targetBox[0] = box->getCoordinate(BoundingBox::XMIN); + targetBox[1] = box->getCoordinate(BoundingBox::XMAX); + targetBox[2] = box->getCoordinate(BoundingBox::YMIN); + targetBox[3] = box->getCoordinate(BoundingBox::YMAX); + targetBox[4] = box->getCoordinate(BoundingBox::ZMIN); + targetBox[5] = box->getCoordinate(BoundingBox::ZMAX); + + std::vector intersectElems; + + tree.getIntersectingElems(targetBox, intersectElems); + + if ( !intersectElems.empty() ) + intersector->intersectCells(targetIdx,intersectElems,result); + } + + delete [] bboxes; + delete [] srcElemIdx; + + // free allocated memory + int ret=intersector->getNumberOfColsOfResMatrix(); + + delete intersector; + + for(unsigned long i = 0 ; i < numSrcElems ; ++i) + { + delete srcElems[i]; + } + for(unsigned long i = 0 ; i < numTargetElems ; ++i) + { + delete targetElems[i]; + } + return ret; + + } +} + +#endif diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am index 1a7fb539e..5c64b6653 100644 --- a/src/INTERP_KERNEL/Makefile.am +++ b/src/INTERP_KERNEL/Makefile.am @@ -51,6 +51,8 @@ Interpolation2D.hxx \ Interpolation2D.txx \ Interpolation3D.hxx \ Interpolation3D.txx \ +Interpolation3D2D.hxx \ +Interpolation3D2D.txx \ Interpolation3DSurf.hxx \ InterpolationOptions.hxx \ InterpolationPlanar.hxx \ @@ -191,6 +193,7 @@ dist_libinterpkernel_la_SOURCES = \ Interpolation2DCurve.cxx \ Interpolation3DSurf.cxx \ Interpolation3D.cxx \ + Interpolation3D2D.cxx \ MeshElement.cxx \ InterpKernelMeshQuality.cxx \ InterpKernelCellSimplify.cxx diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx new file mode 100644 index 000000000..526493e6e --- /dev/null +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.hxx @@ -0,0 +1,65 @@ +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// + +#ifndef __POLYHEDRON3D2DINTERSECTORP0P0_HXX__ +#define __POLYHEDRON3D2DINTERSECTORP0P0_HXX__ + +#include "Intersector3DP0P0.hxx" +#include "SplitterTetra.hxx" +#include "NormalizedUnstructuredMesh.hxx" + +namespace INTERP_KERNEL +{ + + + /** + * \brief Class responsible for calculating intersection between a hexahedron target element and + * the source elements. + * + */ + template + class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0 + { + 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, + SplittingPolicy policy = PLANAR_FACE_5); + + ~Polyhedron3D2DIntersectorP0P0(); + + void intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res); + + private: + void releaseArrays(); + private: + /// pointers to the SplitterTetra objects representing the tetrahedra + /// that result from the splitting of the hexahedron target cell + std::vector< SplitterTetra* > _tetra; + + SplitterTetra2 _split; + + }; +} + +#endif diff --git a/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx new file mode 100644 index 000000000..04d51338a --- /dev/null +++ b/src/INTERP_KERNEL/Polyhedron3D2DIntersectorP0P0.txx @@ -0,0 +1,96 @@ +// Copyright (C) 2007-2011 CEA/DEN, EDF R&D +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +#ifndef __POLYHEDRON3D2DINTERSECTORP0P0_TXX__ +#define __POLYHEDRON3D2DINTERSECTORP0P0_TXX__ + +#include "Polyhedron3D2DIntersectorP0P0.hxx" +#include "Intersector3DP0P0.txx" +#include "MeshUtils.hxx" + +#include "SplitterTetra.txx" + +namespace INTERP_KERNEL +{ + + /** + * Constructor creating object from target cell global number + * The constructor first calculates the necessary nodes, + * (depending on the splitting policy) and then splits the hexahedron into + * tetrahedra, placing these in the internal vector _tetra. + * + * @param targetMesh mesh containing the target elements + * @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) + { + } + + /** + * Destructor. + * Liberates the SplitterTetra objects and potential sub-node points that have been allocated. + * + */ + template + Polyhedron3D2DIntersectorP0P0::~Polyhedron3D2DIntersectorP0P0() + { + releaseArrays(); + } + + template + void Polyhedron3D2DIntersectorP0P0::releaseArrays() + { + for(typename std::vector< SplitterTetra* >::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + delete *iter; + _split.releaseArrays(); + _tetra.clear(); + } + + /** + * Calculates the volume of intersection of an element in the source mesh and the target element + * represented by the object. + * The calculation is performed by calling the corresponding method for + * each SplitterTetra object created by the splitting. + * + * @param targetCell in C mode. + * @param srcCells in C mode. + * + */ + template + void Polyhedron3D2DIntersectorP0P0::intersectCells(ConnType targetCell, const std::vector& srcCells, MyMatrix& res) + { + 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.; + for(typename std::vector*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter) + surface += (*iter)->intersectSourceFace(*iterCellS); + if(surface!=0.) + res[targetCell].insert(std::make_pair(OTT::indFC(*iterCellS), surface)); + } + _split.releaseArrays(); + } +} + +#endif diff --git a/src/INTERP_KERNEL/SplitterTetra.hxx b/src/INTERP_KERNEL/SplitterTetra.hxx index 5127c6d84..3fbe1501d 100644 --- a/src/INTERP_KERNEL/SplitterTetra.hxx +++ b/src/INTERP_KERNEL/SplitterTetra.hxx @@ -173,6 +173,7 @@ namespace INTERP_KERNEL ~SplitterTetra(); double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0); + double intersectSourceFace(typename MyMeshType::MyConnType srcCell); double intersectTetra(const double** tetraCorners); @@ -188,8 +189,10 @@ namespace INTERP_KERNEL // member functions inline void createAffineTransform(const double** corners); inline void checkIsOutside(const double* pt, bool* isOutside) const; + inline void checkIsOutsideSurface(const double* pt, bool* isOutside, const double errTol = DEFAULT_ABS_TOL) const; //dp à supprimer inline void calculateNode(typename MyMeshType::MyConnType globalNodeNum); inline void calculateVolume(TransformedTriangle& tri, const TriangleFaceKey& key); + inline void calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key); /// disallow copying @@ -256,6 +259,21 @@ namespace INTERP_KERNEL isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0); } + //dp à supprimer + template + inline void SplitterTetra::checkIsOutsideSurface(const double* pt, bool* isOutside, const double errTol) const + { + isOutside[0] = isOutside[0] && (pt[0] < -errTol); + isOutside[1] = isOutside[1] && (pt[0] > (1.0 + errTol)); + isOutside[2] = isOutside[2] && (pt[1] < -errTol); + isOutside[3] = isOutside[3] && (pt[1] > (1.0 + errTol)); + isOutside[4] = isOutside[4] && (pt[2] < -errTol); + isOutside[5] = isOutside[5] && (pt[2] > (1.0 + errTol)); + isOutside[6] = isOutside[6] && (1.0 - pt[0] - pt[1] - pt[2] < -errTol); + isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] > (1.0 + errTol)); + } + + /** * Calculates the transformed node with a given global node number. * Gets the coordinates for the node in _src_mesh with the given global number and applies TetraAffineTransform @@ -289,6 +307,22 @@ namespace INTERP_KERNEL _volumes.insert(std::make_pair(key, vol)); } + // TODO DP : adapter les commentaires + // TODO DP : _volume ? + /** + * Calculates the volume contribution from the given TransformedTriangle and stores it with the given key in . + * Calls TransformedTriangle::calculateIntersectionVolume to perform the calculation. + * + * @param tri triangle for which to calculate the volume contribution + * @param key key associated with the face + */ + template + inline void SplitterTetra::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key) + { + const double surf = tri.calculateIntersectionSurface(_t); + _volumes.insert(std::make_pair(key, surf)); + } + template class SplitterTetra2 { diff --git a/src/INTERP_KERNEL/SplitterTetra.txx b/src/INTERP_KERNEL/SplitterTetra.txx index 0372bc4a6..071b7a816 100644 --- a/src/INTERP_KERNEL/SplitterTetra.txx +++ b/src/INTERP_KERNEL/SplitterTetra.txx @@ -306,9 +306,9 @@ namespace INTERP_KERNEL calculateVolume(tri, key1); totalVolume += _volumes[key1]; } else { - // count negative as face has reversed orientation - totalVolume -= _volumes[key1]; - } + // count negative as face has reversed orientation + totalVolume -= _volumes[key1]; + } // local nodes 1, 3, 4 TriangleFaceKey key2 = TriangleFaceKey(faceNodes[0], faceNodes[2], faceNodes[3]); @@ -365,14 +365,257 @@ namespace INTERP_KERNEL { totalVolume = 0.0; } - + LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant()); - // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation + // NB : fault in article, Grandy, [8] : it is the determinant of the inverse transformation // that should be used (which is equivalent to dividing by the determinant) return std::fabs(1.0 / _t->determinant() * totalVolume) ; } + // TODO DP : adapter les commentaires suivants. _volume => _surface ? + /** + * Calculates the volume of intersection of an element in the source mesh and the target element. + * It first calculates the transformation that takes the target tetrahedron into the unit tetrahedron. After that, the + * faces of the source element are triangulated and the calculated transformation is applied + * to each triangle. The algorithm of Grandy, implemented in INTERP_KERNEL::TransformedTriangle is used + * to calculate the contribution to the volume from each triangle. The volume returned is the sum of these contributions + * divided by the determinant of the transformation. + * + * The class will cache the intermediary calculations of transformed nodes of source cells and volumes associated + * with triangulated faces to avoid having to recalculate these. + * + * @param element global number of the source element in C mode. + */ + template + double SplitterTetra::intersectSourceFace(typename MyMeshType::MyConnType element) + { + typedef typename MyMeshType::MyConnType ConnType; + const NumberingPolicy numPol=MyMeshType::My_numPol; + NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT::indFC(element)); + const CellModel& cellModelCell=CellModel::GetCellModel(normCellType); + unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT::indFC(element)) : cellModelCell.getNumberOfNodes(); + + double totalVolume = 0.0; + +#if 1//dp + // calculate the coordinates of the nodes + int *cellNodes=new int[nbOfNodes4Type]; + for(int i = 0;i<(int)nbOfNodes4Type;++i) + { + // 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); + cellNodes[i]=globalNodeNum; + //const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; + } + + // 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 + const double* node = _src_mesh.getCoordinatesPtr()+MyMeshType::MY_SPACEDIM*globalNodeNum; + _nodes[globalNodeNum] = node; + calculateNormalForPolyg(_nodes, nbOfNodes4Type, srcNormal); +#endif + + double faceNormal[3]; + double crossNormals[3]; + for (int iFace = 0; iFace < 4; ++iFace) + { + int decal = iFace * 3; + calculateNormalForTria(_coords + decal, _coords + decal + 1, _coords + decal + 2, 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(_nodes[cellNodes[0]], _nodes[cellNodes[1 + iTri]], _nodes[cellNodes[2 + iTri]], + _coords + decal,_coords + decal + 1,_coords + decal + 2, + 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]); + } + } + break; + } + } +#endif + + //{ could be done on outside? + // check if we have planar tetra element + if(_t->determinant() == 0.0) + { + // tetra is planar + LOG(2, "Planar tetra -- volume 0"); + return 0.0; + } + + // halfspace filtering + bool isOutside[8] = {true, true, true, true, true, true, true, true}; + bool isTargetOutside = false; + + // calculate the coordinates of the nodes +#if 0//dp + int *cellNodes=new int[nbOfNodes4Type]; +#endif + for(int i = 0;i<(int)nbOfNodes4Type;++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); + cellNodes[i]=globalNodeNum; +#else + const int globalNodeNum = cellNodes[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); + } + + checkIsOutsideSurface(_nodes[globalNodeNum], isOutside); + } + + // halfspace filtering check + // NB : might not be beneficial for caching of triangles + for(int i = 0; i < 8; ++i) + { + if(isOutside[i]) + { + isTargetOutside = true; + } + } + + if (!isTargetOutside) + { + // intersect a son with the unit tetra + switch (normCellType) + { + case NORM_TRI3: + { + // create the face key + TriangleFaceKey key = TriangleFaceKey(cellNodes[0], cellNodes[1], cellNodes[2]); + + // calculate the triangle if needed + if (_volumes.find(key) == _volumes.end()) + { + TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1]], _nodes[cellNodes[2]]); + calculateSurface(tri, key); + totalVolume += _volumes[key]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key]; + } + } + break; + + case NORM_QUAD4: //dp pas d'intérêt vis-à-vis du case suivant (à regrouper les deux cases) + + // simple triangulation of faces along a diagonal : + // + // 2 ------ 3 + // | / | + // | / | + // | / | + // | / | + // | / | + // | / | + // 1 ------ 4 + // + //? not sure if this always works + { + // calculate the triangles if needed + + // local nodes 1, 2, 3 + TriangleFaceKey key1 = TriangleFaceKey(cellNodes[0], cellNodes[1], cellNodes[2]); + if (_volumes.find(key1) == _volumes.end()) + { + TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1]], _nodes[cellNodes[2]]); + calculateSurface(tri, key1); + totalVolume += _volumes[key1]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key1]; + } + + // local nodes 1, 3, 4 + TriangleFaceKey key2 = TriangleFaceKey(cellNodes[0], cellNodes[2], cellNodes[3]); + if (_volumes.find(key2) == _volumes.end()) + { + TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[2]], _nodes[cellNodes[3]]); + calculateSurface(tri, key2); + totalVolume += _volumes[key2]; + } + else + { + // count negative as face has reversed orientation + totalVolume -= _volumes[key2]; + } + } + break; + + case NORM_POLYGON: + { + int nbTria = nbOfNodes4Type - 2; // split polygon into nbTria triangles + for (int iTri = 0; iTri < nbTria; ++iTri) + { + TriangleFaceKey key = TriangleFaceKey(cellNodes[0], cellNodes[1 + iTri], cellNodes[2 + iTri]); + if (_volumes.find(key) == _volumes.end()) + { + TransformedTriangle tri(_nodes[cellNodes[0]], _nodes[cellNodes[1 + iTri]], _nodes[cellNodes[2 + iTri]]); + calculateSurface(tri, key); + totalVolume += _volumes[key]; + } + else + { + totalVolume -= _volumes[key]; + } + } + } + break; + + default: + std::cout << "+++ Error : Only elements with triangular and quadratilateral faces are supported at the moment." << std::endl; + assert(false); + } + } + + delete [] cellNodes; + // reset if it is very small to keep the matrix sparse + // is this a good idea? + if(epsilonEqual(totalVolume, 0.0, SPARSE_TRUNCATION_LIMIT)) + { + totalVolume = 0.0; + } + + LOG(2, "Volume = " << totalVolume << ", det= " << _t->determinant()); + + return totalVolume; + } + /** * Calculates the volume of intersection of this tetrahedron with another one. */ @@ -643,7 +886,12 @@ namespace INTERP_KERNEL int conn[4]; for(int j = 0; j < 4; ++j) { - nodes[j] = getCoordsOfSubNode2(subZone[ SPLIT_NODES_6[4*i+j] ],conn[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); diff --git a/src/INTERP_KERNEL/TransformedTriangle.cxx b/src/INTERP_KERNEL/TransformedTriangle.cxx index 661d71e65..538b13a18 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.cxx +++ b/src/INTERP_KERNEL/TransformedTriangle.cxx @@ -19,7 +19,9 @@ #include "TransformedTriangle.hxx" #include "VectorUtils.hxx" - +#if 1//dp +#include "TetraAffineTransform.hxx" +#endif #include #include #include @@ -142,7 +144,7 @@ namespace INTERP_KERNEL * Destructor * * Deallocates the memory used to store the points of the polygons. - * This memory is allocated in calculateIntersectionPolygons(). + * This memory is allocated in calculateIntersectionAndProjectionPolygons(). */ TransformedTriangle::~TransformedTriangle() { @@ -197,7 +199,7 @@ namespace INTERP_KERNEL //sign = sign > 0.0 ? 1.0 : -1.0; LOG(2, "-- Calculating intersection polygons ... "); - calculateIntersectionPolygons(); + calculateIntersectionAndProjectionPolygons(); double barycenter[3]; @@ -230,6 +232,40 @@ namespace INTERP_KERNEL } + /** + * Calculates the volume of intersection between the triangle and the + * unit tetrahedron. + * + * @return volume of intersection of this triangle with unit tetrahedron, + * as described in Grandy + * + */ + double TransformedTriangle::calculateIntersectionSurface(TetraAffineTransform* tat) + { + // check first that we are not below z - plane + if(isTriangleBelowTetraeder()) + { + LOG(2, " --- Triangle is below tetraeder - V = 0.0"); + return 0.0; + } + + LOG(2, "-- Calculating intersection polygon ... "); + calculateIntersectionPolygon(); + + _volume = 0.; + if(_polygonA.size() > 2) { + double barycenter[3]; + calculatePolygonBarycenter(A, barycenter); + sortIntersectionPolygon(A, barycenter); + const int nbPoints = _polygonA.size(); + for(int i = 0 ; i < nbPoints ; ++i) + tat->reverseApply(_polygonA[i], _polygonA[i]); + _volume = calculateSurfacePolygon(); + } + + return _volume; + } + // ---------------------------------------------------------------------------------- // TransformedTriangle PRIVATE // ---------------------------------------------------------------------------------- @@ -242,7 +278,7 @@ namespace INTERP_KERNEL * intersection polygon B. * */ - void TransformedTriangle::calculateIntersectionPolygons() + void TransformedTriangle::calculateIntersectionAndProjectionPolygons() { assert(_polygonA.size() == 0); assert(_polygonB.size() == 0); @@ -437,10 +473,137 @@ namespace INTERP_KERNEL } + /** + * Calculates the intersection polygon A, performing the intersection tests + * and storing the corresponding point in the vector _polygonA. + * + * @post _polygonA contains the intersection polygon A. + * + */ + void TransformedTriangle::calculateIntersectionPolygon() + { + assert(_polygonA.size() == 0); + // avoid reallocations in push_back() by pre-allocating enough memory + // we should never have more than 20 points + _polygonA.reserve(20); + // -- surface intersections + // surface - edge + for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { + if(testSurfaceEdgeIntersection(edge)) + { + double* ptA = new double[3]; + calcIntersectionPtSurfaceEdge(edge, ptA); + _polygonA.push_back(ptA); + LOG(3,"Surface-edge : " << vToStr(ptA) << " added to A "); + } + } + + // -- segment intersections + for(TriSegment seg = PQ ; seg < NO_TRI_SEGMENT ; seg = TriSegment(seg + 1)) + { + + bool isZero[NO_DP]; + + // check beforehand which double-products are zero + for(DoubleProduct dp = C_YZ; dp < NO_DP; dp = DoubleProduct(dp + 1)) + { + isZero[dp] = (calcStableC(seg, dp) == 0.0); + } + + // segment - facet + for(TetraFacet facet = OYZ ; facet < NO_TET_FACET ; facet = TetraFacet(facet + 1)) + { + // is this test worth it? + const bool doTest = + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet]] && + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 1]] && + !isZero[DP_FOR_SEG_FACET_INTERSECTION[3*facet + 2]]; + + if(doTest && testSegmentFacetIntersection(seg, facet)) + { + double* ptA = new double[3]; + calcIntersectionPtSegmentFacet(seg, facet, ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-facet : " << vToStr(ptA) << " added to A"); + } + } + + // segment - edge + for(TetraEdge edge = OX ; edge <= ZX ; edge = TetraEdge(edge + 1)) + { + const DoubleProduct edge_dp = DoubleProduct(edge); + + if(isZero[edge_dp] && testSegmentEdgeIntersection(seg, edge)) + { + double* ptA = new double[3]; + calcIntersectionPtSegmentEdge(seg, edge, ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-edge : " << vToStr(ptA) << " added to A"); + } + } + + // segment - corner + for(TetraCorner corner = O ; corner < NO_TET_CORNER ; corner = TetraCorner(corner + 1)) + { + const bool doTest = + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner] )] && + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+1] )] && + isZero[DoubleProduct( EDGES_FOR_CORNER[3*corner+2] )]; + + if(doTest && testSegmentCornerIntersection(seg, corner)) + { + double* ptA = new double[3]; + copyVector3(&COORDS_TET_CORNER[3 * corner], ptA); + _polygonA.push_back(ptA); + LOG(3,"Segment-corner : " << vToStr(ptA) << " added to A"); + } + } + + } + + // inclusion tests + for(TriCorner corner = P ; corner < NO_TRI_CORNER ; corner = TriCorner(corner + 1)) + { + // { XYZ - inclusion only possible if in Tetrahedron? + // tetrahedron + if(testCornerInTetrahedron(corner)) + { + double* ptA = new double[3]; + copyVector3(&_coords[5*corner], ptA); + _polygonA.push_back(ptA); + LOG(3,"Inclusion tetrahedron : " << vToStr(ptA) << " added to A"); + } + + } + + } + + //TODO DP : commentaires + double TransformedTriangle::calculateSurfacePolygon() + { + const int nbPoints = _polygonA.size(); + double pdt[3]; + double sum[3] = {0., 0., 0.}; + + 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) + + cross(ptCurr, ptNext, pdt); + add(pdt, sum); + } + + const double surface = norm(sum) * 0.5; + LOG(2,"Surface is " << surface); + return surface; + } + /** * Calculates the barycenters of the given intersection polygon. * - * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() + * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons() * * @param poly one of the two intersection polygons * @param barycenter array of three doubles where barycenter is stored @@ -477,7 +640,7 @@ namespace INTERP_KERNEL /** * Sorts the given intersection polygon in circular order around its barycenter. - * @pre the intersection polygons have been calculated with calculateIntersectionPolygons() + * @pre the intersection polygons have been calculated with calculateIntersectionAndProjectionPolygons() * @post the vertices in _polygonA and _polygonB are sorted in circular order around their * respective barycenters * @@ -537,7 +700,7 @@ namespace INTERP_KERNEL /** * Calculates the volume between the given polygon and the z = 0 plane. * - * @pre the intersection polygones have been calculated with calculateIntersectionPolygons(), + * @pre the intersection polygones have been calculated with calculateIntersectionAndProjectionPolygons(), * and they have been sorted in circular order with sortIntersectionPolygons(void) * * @param poly one of the two intersection polygons diff --git a/src/INTERP_KERNEL/TransformedTriangle.hxx b/src/INTERP_KERNEL/TransformedTriangle.hxx index 84877b9b9..2dba44370 100644 --- a/src/INTERP_KERNEL/TransformedTriangle.hxx +++ b/src/INTERP_KERNEL/TransformedTriangle.hxx @@ -46,6 +46,9 @@ 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 @@ -137,6 +140,7 @@ namespace INTERP_KERNEL ~TransformedTriangle(); double calculateIntersectionVolume(); + double calculateIntersectionSurface(TetraAffineTransform* tat); void dumpCoords() const; @@ -155,7 +159,7 @@ namespace INTERP_KERNEL // ---------------------------------------------------------------------------------- // High-level methods called directly by calculateIntersectionVolume() // ---------------------------------------------------------------------------------- - void calculateIntersectionPolygons(); + void calculateIntersectionAndProjectionPolygons(); void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter); @@ -163,6 +167,13 @@ namespace INTERP_KERNEL double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter); + // ---------------------------------------------------------------------------------- + // High-level methods called directly by calculateIntersectionSurface() + // ---------------------------------------------------------------------------------- + void calculateIntersectionPolygon(); + + double calculateSurfacePolygon(); + // ---------------------------------------------------------------------------------- // Detection of degenerate triangles // ---------------------------------------------------------------------------------- diff --git a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx index 6aecef649..03737720a 100644 --- a/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx +++ b/src/INTERP_KERNEL/UnitTetraIntersectionBary.cxx @@ -94,7 +94,7 @@ namespace INTERP_KERNEL *((TransformedTriangle*)this) = triangle; // copy triangle fields _polygonA.clear(); _polygonB.clear(); - calculateIntersectionPolygons(); + calculateIntersectionAndProjectionPolygons(); if (this->_polygonA.size() < 3) return; calculatePolygonBarycenter(A, _barycenterA); diff --git a/src/INTERP_KERNEL/VectorUtils.hxx b/src/INTERP_KERNEL/VectorUtils.hxx index bc35d333c..1fcf93e0e 100644 --- a/src/INTERP_KERNEL/VectorUtils.hxx +++ b/src/INTERP_KERNEL/VectorUtils.hxx @@ -80,6 +80,35 @@ namespace INTERP_KERNEL return ss.str(); } + + //dp : à supprimer + /** + * Subtracts two a double[3] - vectors and store the result in res + * + * @param v1 vector v1 + * @param v2 vector v2 + * @param res vector in which to store the result v1 - v2. + */ + inline void subtract(const double* v1, const double* v2, double* res) + { + res[0] = v1[0] - v2[0]; + res[1] = v1[1] - v2[1]; + res[2] = v1[2] - v2[2]; + } + + /** + * Adds a double[3] - vector to another one. + * + * @param v vector v + * @param res vector in which to store the result res + v. + */ + inline void add(const double* v, double* res) + { + res[0] += v[0]; + res[1] += v[1]; + res[2] += v[2]; + } + /** * Calculates the cross product of two double[3] - vectors. * diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx index 6fbe69a54..5b4025418 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx @@ -36,6 +36,7 @@ namespace ParaMEDMEM { CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest); //MEDCouplingBasicsTest1.cxx +#if 0//dp CPPUNIT_TEST( testArray ); CPPUNIT_TEST( testArray2 ); CPPUNIT_TEST( testArray3 ); @@ -252,8 +253,10 @@ 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 ); @@ -274,13 +277,18 @@ namespace ParaMEDMEM CPPUNIT_TEST( test3DSurfInterpP1P0Bary_1 ); CPPUNIT_TEST( test3DInterpP1P0Bary_1 ); CPPUNIT_TEST( test3DTo1DInterpP0P0PL_1 ); +#endif + + CPPUNIT_TEST( test3D2DInterpP0P0_1 ); +#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 @@ -531,6 +539,10 @@ namespace ParaMEDMEM void test2DCurveInterpP1P0_1(); void test2DCurveInterpP1P1_1(); +#if 1//dp + void test3D2DInterpP0P0_1(); +#endif + public: static MEDCouplingUMesh *build3DSourceMesh_2(); static MEDCouplingUMesh *build3DTargetMesh_2(); @@ -570,6 +582,12 @@ namespace ParaMEDMEM static MEDCouplingUMesh *build3DMultiTypes_1(); static MEDCouplingUMesh *buildHexa8Mesh_1(); static MEDCouplingUMesh *buildPointe_1(MEDCouplingUMesh *&m1); + +#if 1//dp + static MEDCouplingUMesh *build3D2DSourceMesh_1(); + static MEDCouplingUMesh *build3D2DTargetMesh_1(); +#endif + static DataArrayDouble *buildCoordsForMultiTypes_1(); static MEDCouplingMultiFields *buildMultiFields_1(); static std::vector buildMultiFields_2(); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx index 936adb68f..456343d41 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest0.cxx @@ -1173,3 +1173,85 @@ double MEDCouplingBasicsTest::sumAll(const std::vector< std::map >& ret+=(*iter2).second; return ret; } + +#if 1//dp + +#if 1 +MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1() +{ + double sourceCoords[42]={-12., 6., 10., -12.,10., 6., -16.,10., 10., + -20., 0., 0., -12., 0., 0., -12., 0., -4., -20.,0.,-4., + -20., 0., 10., -12., 0., 10., -20.,10., 10., + -25., 5., -5., 5., 5., -5., 5., 5., 25., -25.,5.,25.}; + int sourceConn[14]={0,1,2, 3,4,5,6, 7,8,9, 10,11,12,13}; + MEDCouplingUMesh *sourceMesh=MEDCouplingUMesh::New(); + sourceMesh->setMeshDimension(2); + sourceMesh->allocateCells(1); + //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn); + //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+3); + //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3 ,3,sourceConn+7); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,sourceConn+10); + sourceMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(14,3); + std::copy(sourceCoords,sourceCoords+42,myCoords->getPointer()); + sourceMesh->setCoords(myCoords); + myCoords->decrRef(); + return sourceMesh; +} + +MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh_1() +{ + 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(1); + //sourceMesh->insertNextCell(INTERP_KERNEL::NORM_TETRA4,4,sourceConn); + sourceMesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8,8,sourceConn + 4); + sourceMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(9,3); + std::copy(sourceCoords,sourceCoords+27,myCoords->getPointer()); + sourceMesh->setCoords(myCoords); + myCoords->decrRef(); + return sourceMesh; +} + +#else // détruire dessous + +MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DTargetMesh_1() +{ + 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); + targetMesh->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(4,3); + std::copy(targetCoords,targetCoords+12,myCoords->getPointer()); + targetMesh->setCoords(myCoords); + myCoords->decrRef(); + return targetMesh; +} + +MEDCouplingUMesh *MEDCouplingBasicsTest::build3D2DSourceMesh_1() +{ + 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->finishInsertingCells(); + DataArrayDouble *myCoords=DataArrayDouble::New(); + myCoords->alloc(3,3); + std::copy(targetCoords,targetCoords+9,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 03984bc64..02d97c2d1 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTestInterp.cxx @@ -25,6 +25,7 @@ #include "Interpolation2D.txx" #include "Interpolation3DSurf.hxx" #include "Interpolation3D.txx" +#include "Interpolation3D2D.txx" #include "InterpolationCC.txx" #include "InterpolationCU.txx" #include "Interpolation2DCurve.hxx" @@ -2265,3 +2266,69 @@ void MEDCouplingBasicsTest::test2DCurveInterpP1P1_1() sourceMesh->decrRef(); targetMesh->decrRef(); } + +#if 1//dp +#if 0 +void MEDCouplingBasicsTest::test3D2DInterpP0P0_1() +{ + MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3D2DTargetMesh_1(); + + 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::PLANAR_FACE_5, INTERP_KERNEL::PLANAR_FACE_6, INTERP_KERNEL::GENERAL_24, INTERP_KERNEL::GENERAL_48 }; + for ( int i = 0; i < 4; ++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); + } + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} +#else +void MEDCouplingBasicsTest::test3D2DInterpP0P0_1() +{ + MEDCouplingUMesh *sourceMesh=build3D2DSourceMesh_1(); + MEDCouplingUMesh *targetMesh=build3D2DTargetMesh_1(); + + 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 }; + //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 ) + { + myInterpolator.setSplittingPolicy( sp[i] ); + res.clear(); + myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,res,"P0P0"); + + CPPUNIT_ASSERT_EQUAL(1,(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(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[0][0],1e-12); + + //CPPUNIT_ASSERT_DOUBLES_EQUAL(8.e6,sumAll(res),1e-7); + } + //clean up + sourceMesh->decrRef(); + targetMesh->decrRef(); +} +#endif +#endif