--- /dev/null
+// 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<Interpolation3D2D>(io)
+ {
+ }
+}
--- /dev/null
+// 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<Interpolation3D2D>
+ {
+ public:
+ Interpolation3D2D();
+ Interpolation3D2D(const InterpolationOptions& io);
+ template<class MyMeshType, class MatrixType>
+ int interpolateMeshes(const MyMeshType& srcMesh, const MyMeshType& targetMesh, MatrixType& result, const char *method);
+ private:
+ SplittingPolicy _splitting_policy;
+ };
+}
+
+#endif
--- /dev/null
+// 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 <stack>
+
+#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<class MyMeshType, class MatrixType>
+ 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<MeshElement<ConnType>*> srcElems(numSrcElems);
+ std::vector<MeshElement<ConnType>*> targetElems(numTargetElems);
+
+ std::map<MeshElement<ConnType>*, int> indices;
+
+ for(unsigned long i = 0 ; i < numSrcElems ; ++i)
+ srcElems[i] = new MeshElement<ConnType>(i, srcMesh);
+
+ for(unsigned long i = 0 ; i < numTargetElems ; ++i)
+ targetElems[i] = new MeshElement<ConnType>(i, targetMesh);
+
+ Intersector3D<MyMeshType,MatrixType>* intersector=0;
+ std::string methC = InterpolationOptions::filterInterpolationMethod(method);
+ if(methC=="P0P0")
+ {
+ switch(InterpolationOptions::getIntersectionType())
+ {
+ case Triangulation:
+ intersector=new Polyhedron3D2DIntersectorP0P0<MyMeshType,MatrixType>(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<ConnType> 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
Interpolation2D.txx \
Interpolation3D.hxx \
Interpolation3D.txx \
+Interpolation3D2D.hxx \
+Interpolation3D2D.txx \
Interpolation3DSurf.hxx \
InterpolationOptions.hxx \
InterpolationPlanar.hxx \
Interpolation2DCurve.cxx \
Interpolation3DSurf.cxx \
Interpolation3D.cxx \
+ Interpolation3D2D.cxx \
MeshElement.cxx \
InterpKernelMeshQuality.cxx \
InterpKernelCellSimplify.cxx
--- /dev/null
+// 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 MyMeshType, class MyMatrix>
+ class Polyhedron3D2DIntersectorP0P0 : public Intersector3DP0P0<MyMeshType,MyMatrix>
+ {
+ 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<ConnType>& 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<MyMeshType>* > _tetra;
+
+ SplitterTetra2<MyMeshType> _split;
+
+ };
+}
+
+#endif
--- /dev/null
+// 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<class MyMeshType, class MyMatrix>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::Polyhedron3D2DIntersectorP0P0(const MyMeshType& targetMesh, const MyMeshType& srcMesh,
+ SplittingPolicy policy)
+ : Intersector3DP0P0<MyMeshType,MyMatrix>(targetMesh,srcMesh),_split(targetMesh,srcMesh,policy)
+ {
+ }
+
+ /**
+ * Destructor.
+ * Liberates the SplitterTetra objects and potential sub-node points that have been allocated.
+ *
+ */
+ template<class MyMeshType, class MyMatrix>
+ Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::~Polyhedron3D2DIntersectorP0P0()
+ {
+ releaseArrays();
+ }
+
+ template<class MyMeshType, class MyMatrix>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::releaseArrays()
+ {
+ for(typename std::vector< SplitterTetra<MyMeshType>* >::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<class MyMeshType, class MyMatrix>
+ void Polyhedron3D2DIntersectorP0P0<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
+ {
+ int nbOfNodesT=Intersector3D<MyMeshType,MyMatrix>::_target_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::indFC(targetCell));
+ releaseArrays();
+ _split.splitTargetCell(targetCell,nbOfNodesT,_tetra);
+ for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
+ {
+ double surface = 0.;
+ for(typename std::vector<SplitterTetra<MyMeshType>*>::iterator iter = _tetra.begin(); iter != _tetra.end(); ++iter)
+ surface += (*iter)->intersectSourceFace(*iterCellS);
+ if(surface!=0.)
+ res[targetCell].insert(std::make_pair(OTT<ConnType,numPol>::indFC(*iterCellS), surface));
+ }
+ _split.releaseArrays();
+ }
+}
+
+#endif
~SplitterTetra();
double intersectSourceCell(typename MyMeshType::MyConnType srcCell, double* baryCentre=0);
+ double intersectSourceFace(typename MyMeshType::MyConnType srcCell);
double intersectTetra(const double** tetraCorners);
// 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
isOutside[7] = isOutside[7] && (1.0 - pt[0] - pt[1] - pt[2] >= 1.0);
}
+ //dp à supprimer
+ template<class MyMeshType>
+ inline void SplitterTetra<MyMeshType>::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
_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<class MyMeshType>
+ inline void SplitterTetra<MyMeshType>::calculateSurface(TransformedTriangle& tri, const TriangleFaceKey& key)
+ {
+ const double surf = tri.calculateIntersectionSurface(_t);
+ _volumes.insert(std::make_pair(key, surf));
+ }
+
template<class MyMeshTypeT, class MyMeshTypeS=MyMeshTypeT>
class SplitterTetra2
{
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]);
{
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<class MyMeshType>
+ double SplitterTetra<MyMeshType>::intersectSourceFace(typename MyMeshType::MyConnType element)
+ {
+ typedef typename MyMeshType::MyConnType ConnType;
+ const NumberingPolicy numPol=MyMeshType::My_numPol;
+ NormalizedCellType normCellType=_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(element));
+ const CellModel& cellModelCell=CellModel::GetCellModel(normCellType);
+ unsigned nbOfNodes4Type=cellModelCell.isDynamic() ? _src_mesh.getNumberOfNodesOfElement(OTT<ConnType,numPol>::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<ConnType,numPol>::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<double> 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<MyMeshType,MyMatrix>::_dim_caracteristic,
+ DEFAULT_ABS_TOL); //dp PlanarIntersector<MyMeshType,MyMatrix>::_precision);
+ ConnType nb_inter=((ConnType)inter.size())/2;
+ if(nb_inter >3) inter=reconstruct_polygon(inter);
+ for(ConnType i = 1; i<nb_inter-1; i++)
+ {
+ INTERP_KERNEL::crossprod<2>(&inter[0],&inter[2*i],&inter[2*(i+1)],area);
+ totalVolume +=0.5*fabs(area[0]);
+ }
+ }
+ 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<ConnType,numPol>::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.
*/
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<MyMeshTypeS>* t = new SplitterTetra<MyMeshTypeS>(_src_mesh, nodes,conn);
tetra.push_back(t);
#include "TransformedTriangle.hxx"
#include "VectorUtils.hxx"
-
+#if 1//dp
+#include "TetraAffineTransform.hxx"
+#endif
#include <iostream>
#include <fstream>
#include <cassert>
* 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()
{
//sign = sign > 0.0 ? 1.0 : -1.0;
LOG(2, "-- Calculating intersection polygons ... ");
- calculateIntersectionPolygons();
+ calculateIntersectionAndProjectionPolygons();
double barycenter[3];
}
+ /**
+ * 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
// ----------------------------------------------------------------------------------
* intersection polygon B.
*
*/
- void TransformedTriangle::calculateIntersectionPolygons()
+ void TransformedTriangle::calculateIntersectionAndProjectionPolygons()
{
assert(_polygonA.size() == 0);
assert(_polygonB.size() == 0);
}
+ /**
+ * 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
/**
* 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
*
/**
* 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
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
~TransformedTriangle();
double calculateIntersectionVolume();
+ double calculateIntersectionSurface(TetraAffineTransform* tat);
void dumpCoords() const;
// ----------------------------------------------------------------------------------
// High-level methods called directly by calculateIntersectionVolume()
// ----------------------------------------------------------------------------------
- void calculateIntersectionPolygons();
+ void calculateIntersectionAndProjectionPolygons();
void calculatePolygonBarycenter(const IntersectionPolygon poly, double* barycenter);
double calculateVolumeUnderPolygon(IntersectionPolygon poly, const double* barycenter);
+ // ----------------------------------------------------------------------------------
+ // High-level methods called directly by calculateIntersectionSurface()
+ // ----------------------------------------------------------------------------------
+ void calculateIntersectionPolygon();
+
+ double calculateSurfacePolygon();
+
// ----------------------------------------------------------------------------------
// Detection of degenerate triangles
// ----------------------------------------------------------------------------------
*((TransformedTriangle*)this) = triangle; // copy triangle fields
_polygonA.clear();
_polygonB.clear();
- calculateIntersectionPolygons();
+ calculateIntersectionAndProjectionPolygons();
if (this->_polygonA.size() < 3)
return;
calculatePolygonBarycenter(A, _barycenterA);
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.
*
{
CPPUNIT_TEST_SUITE(MEDCouplingBasicsTest);
//MEDCouplingBasicsTest1.cxx
+#if 0//dp
CPPUNIT_TEST( testArray );
CPPUNIT_TEST( testArray2 );
CPPUNIT_TEST( testArray3 );
CPPUNIT_TEST( testInterpolationCU1D );
CPPUNIT_TEST( testInterpolationCU2D );
CPPUNIT_TEST( testInterpolationCU3D );
+#endif
CPPUNIT_TEST( test3DInterpP0P0_1 );
+#if 0 //dp
CPPUNIT_TEST( test3DInterpP0P0PL_1 );
CPPUNIT_TEST( test3DInterpP0P0PL_2 );
CPPUNIT_TEST( test3DInterpP0P0PL_3 );
CPPUNIT_TEST( test3DSurfInterpP1P0Bary_1 );
CPPUNIT_TEST( test3DInterpP1P0Bary_1 );
CPPUNIT_TEST( test3DTo1DInterpP0P0PL_1 );
+#endif
+
+ CPPUNIT_TEST( test3D2DInterpP0P0_1 );
+#if 0//dp
CPPUNIT_TEST( test1DInterp_1 );
CPPUNIT_TEST( test2DCurveInterpP0P0_1 );
CPPUNIT_TEST( test2DCurveInterpP0P0_2 );
CPPUNIT_TEST( test2DCurveInterpP0P1_1 );
CPPUNIT_TEST( test2DCurveInterpP1P0_1 );
CPPUNIT_TEST( test2DCurveInterpP1P1_1 );
+#endif
CPPUNIT_TEST_SUITE_END();
public:
//MEDCouplingBasicsTest1.cxx
void test2DCurveInterpP1P0_1();
void test2DCurveInterpP1P1_1();
+#if 1//dp
+ void test3D2DInterpP0P0_1();
+#endif
+
public:
static MEDCouplingUMesh *build3DSourceMesh_2();
static MEDCouplingUMesh *build3DTargetMesh_2();
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<MEDCouplingFieldDouble *> buildMultiFields_2();
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
#include "Interpolation2D.txx"
#include "Interpolation3DSurf.hxx"
#include "Interpolation3D.txx"
+#include "Interpolation3D2D.txx"
#include "InterpolationCC.txx"
#include "InterpolationCU.txx"
#include "Interpolation2DCurve.hxx"
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<std::map<int,double> > 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<std::map<int,double> > 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
