Remapper: added new P1P1 method: MappedBarycentric. See doc

for more infos.

diff --git a/doc/user/doxygen/doxfiles/reference/interpolation/intersectors.dox b/doc/user/doxygen/doxfiles/reference/interpolation/intersectors.dox
index 1729060bd3609d0b5cd00840a4e8459c6d7c4115..8d88abda25366da4c5f1b2bcaf9c2ab26884c43c 100644 (file)
--- a/doc/user/doxygen/doxfiles/reference/interpolation/intersectors.dox
+++ b/doc/user/doxygen/doxfiles/reference/interpolation/intersectors.dox
@@ -11,7 +11,8 @@ Before reading on, remember the definition of a \ref glossary "P0 and P1 field".
 
 - \subpage intersec-specifics
 - \subpage interpkernelGeo2D
-- \subpage barycoords (used in some P1 intersectors/locators) 
+- \subpage barycoords (used in some P1 intersectors/locators)
+- \subpage mapped_bary (used in some P1P1 intersectors) 
 
 Some implementation details of the C++ code can also be found here: \ref interpkernel 
 

diff --git a/doc/user/doxygen/doxfiles/reference/interpolation/mapped_bary.dox b/doc/user/doxygen/doxfiles/reference/interpolation/mapped_bary.dox
new file mode 100644 (file)
index 0000000..c724d85
--- /dev/null
+++ b/doc/user/doxygen/doxfiles/reference/interpolation/mapped_bary.dox
@@ -0,0 +1,61 @@
+/*!
+\page mapped_bary Mapped barycentric coordinates algorithm
+
+Mapped barycentric intersection type ('MappedBarycentric') can be selected in space dim 2 (resp. 3) when
+working with quadrangle only (resp. hexaedrons only).
+
+It can only be used for P1P1 projection: for any point P within the quadrangle or the hexaedron, the 
+set of reduced coordinates is computed (x, y, z all comprised between 0 and 1).
+Then the field value at P is computed using the usual form functions of finite element method 
+((1-x)*(1-y), x*(1-y), (1-x)*y and x*y in 2D for example). 
+
+The algorithm used to compute the reduced coordinates differs in dim 2 and dim 3.
+
+\section mapped_bary2d Dimension 2
+
+Let O, A, B, C the four points of the quadrangle, clockwise. Without loss of generality 
+O is assumed to be the origin.
+A point P within the quadrangle is identified with vector OP and simply denoted P.
+
+A suitable mapping is such that, if (x,y) is the couple of reduced coordinates (with x and y both in [0,1])
+ of a point P, then:
+\f[   \mathbf{P} = x\mathbf{C} + y\mathbf{A} + xy(\mathbf{B}- \mathbf{A}-\mathbf{C}) \f]
+
+This forms is the simplest one having a gradient which x component is constant in x, and similarly in y. 
+Furthermore the reduced coordinates 
+(0,0) (resp. (0,1), (1,0), and (1,1)) map to point O (resp. A, B, and C).
+
+Calling \f$\mathbf{N} = \mathbf{B}-\mathbf{A}-\mathbf{C}\f$ and developping for the 2 compos:
+
+\f[   p_x = C_x x + A_x y + N_x xy \f]
+\f[   p_y = C_y x + A_y y + N_y xy \f]
+
+Solving the first eq for x:
+\f[ x = \frac{p_x - y A_x }{C_x+yN_x}  \f]
+
+and injecting in second eq:
+\f[   (A_yN_x -N_yA_x)y^2 + (-p_yN_x -A_xC_y +A_yC_x+N_yp_x)y + (p_x C_y-p_yC_x)=0   \f] 
+solved in:
+\f[  y = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2a} \f]
+where a, b and c are the coeffs of the 2nd order eq.
+
+This gives two possible couples of solution among which only one is in \f$[0,1]^2\f$.
+
+In case where \f$A_yN_x -N_yA_x = 0 \f$ we have a degenerated unique solution for $y$
+\f[ y = \frac{c}{b} \f]
+
+\subsection{Rectangle}
+
+Finally it is worth puting aside the case \f$\mathbf{N} = 0\f$ (rectangle), which boils down to solving an ordinary
+2-unknows system:
+\f[ x = \frac{p_x A_y - p_y A_x}{C_x A_y - C_y A_x}, y = \frac{C_x p_y-C_y p_x}{C_x A_y - C_y A_x} \f]
+
+
+\section mapped_bary3d Dimension 3  
+
+In three dimensions, adopting the same approach as above would lead to a 4th order equation to solve.
+A simpler approach has been chosen: the distance to each pair of parallel faces in the hexaedron is computed.   
+The ratios to the sum of the two distances is computed giving again a number between 0 and 1 for each of 
+the 3 directions.
+
+*/

diff --git a/src/INTERP_KERNEL/InterpKernelMatrix.hxx b/src/INTERP_KERNEL/InterpKernelMatrix.hxx
old mode 100755 (executable)
new mode 100644 (file)

diff --git a/src/INTERP_KERNEL/Interpolation1D.hxx b/src/INTERP_KERNEL/Interpolation1D.hxx
old mode 100755 (executable)
new mode 100644 (file)

diff --git a/src/INTERP_KERNEL/Interpolation2D.hxx b/src/INTERP_KERNEL/Interpolation2D.hxx
old mode 100755 (executable)
new mode 100644 (file)

diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx
index eecea5607bc2faeb54d421bd07325d3535d3190c..c78f2a059b5f2e412b20579b15032134c4d145e6 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation3D.txx
+++ b/src/INTERP_KERNEL/Interpolation3D.txx
@@ -35,6 +35,7 @@
 #include "PolyhedronIntersectorP1P1.txx"
 #include "PointLocator3DIntersectorP1P1.txx"
 #include "Barycentric3DIntersectorP1P1.txx"
+#include "MappedBarycentric3DIntersectorP1P1.txx"
 #include "Log.hxx"
 // If defined, use recursion to traverse the binary search tree, else use the BBTree class
 //#define USE_RECURSIVE_BBOX_FILTER
@@ -154,8 +155,11 @@ namespace INTERP_KERNEL
           case Barycentric:
             intersector=new Barycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
             break;
+          case MappedBarycentric:
+            intersector=new MappedBarycentric3DIntersectorP1P1<MyMeshType,MatrixType>(targetMesh, srcMesh, getPrecision());
+            break;
           default:
-            throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle or PointLocator.");
+            throw INTERP_KERNEL::Exception("Invalid 3D intersection type for P1P1 interp specified : must be Triangle, PointLocator, Barycentric or MappedBarycentric.");
           }
       }
     else

diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx
index 74fe301c943634499764e1174d8973061cfce987..4187a69d5c7e671642cf6f39aa4cd0ac4cc1b56c 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationOptions.hxx
+++ b/src/INTERP_KERNEL/InterpolationOptions.hxx
@@ -28,7 +28,7 @@
 
 namespace INTERP_KERNEL
 {
-  typedef enum { Triangulation, Convex, Geometric2D, PointLocator, Barycentric, BarycentricGeo2D } IntersectionType;
+  typedef enum { Triangulation, Convex, Geometric2D, PointLocator, Barycentric, BarycentricGeo2D, MappedBarycentric } IntersectionType;
   
   /*!
    * Class defining the options for all interpolation algorithms used in the \ref remapper "remapper" and

diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx
old mode 100755 (executable)
new mode 100644 (file)

diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx
index b6ff71effe0040a363c879a98bd2f77a74c42432..0dc6e42b886cd5cf98b87a33bc92bca7d3b9037d 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationPlanar.txx
+++ b/src/INTERP_KERNEL/InterpolationPlanar.txx
@@ -40,6 +40,8 @@
 #include "PlanarIntersectorP1P0PL.txx"
 #include "PlanarIntersectorP1P1PL.hxx"
 #include "PlanarIntersectorP1P1PL.txx"
+#include "MappedBarycentric2DIntersectorP1P1.hxx"
+#include "MappedBarycentric2DIntersectorP1P1.txx"
 #include "VectorUtils.hxx"
 #include "BBTree.txx"
 
@@ -358,8 +360,16 @@ namespace INTERP_KERNEL
                                                                            InterpolationOptions::getPrecision(),
                                                                            InterpolationOptions::getOrientation());
             break;
+          case MappedBarycentric:
+            intersector=new MappedBarycentric2DIntersectorP1P1<MyMeshType,MatrixType>(myMeshT, myMeshS, _dim_caracteristic,
+                                                                                     InterpolationOptions::getMaxDistance3DSurfIntersect(),
+                                                                                     InterpolationOptions::getMinDotBtwPlane3DSurfIntersect(),
+                                                                                     InterpolationOptions::getMedianPlane(),
+                                                                                     InterpolationOptions::getPrecision(),
+                                                                                     InterpolationOptions::getOrientation());
+            break;
           default:
-            throw INTERP_KERNEL::Exception("For P1P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator !");
+            throw INTERP_KERNEL::Exception("For P1P1 planar interpolation possibities are : Triangulation, Convex, Geometric2D, PointLocator, MappedBarycentric !");
           }
       }
     else

diff --git a/src/INTERP_KERNEL/InterpolationUtils.hxx b/src/INTERP_KERNEL/InterpolationUtils.hxx
index 8f7eaffa6ac16125248251c86cc3379873a6db3e..f8209fd68605d548844bc7647c75866592a6d716 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationUtils.hxx
+++ b/src/INTERP_KERNEL/InterpolationUtils.hxx
@@ -23,6 +23,7 @@
 
 #include "INTERPKERNELDefines.hxx"
 #include "InterpKernelException.hxx"
+#include "VolSurfUser.hxx"
 
 #include "NormalizedUnstructuredMesh.hxx"
 
@@ -413,7 +414,7 @@ namespace INTERP_KERNEL
   }
 
   /*!
-   * Calculate barycentric coordinates of a point p with respect to triangle or tetra verices.
+   * Calculate barycentric coordinates of a point p with respect to triangle or tetra vertices.
    * This method makes 2 assumptions :
    *    - this is a simplex
    *    - spacedim == meshdim. For TRI3 and TRI6 spaceDim is expected to be equal to 2 and for TETRA4 spaceDim is expected to be equal to 3.
@@ -535,6 +536,136 @@ namespace INTERP_KERNEL
       }
   }
 
+  /*!
+   * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices.
+   * This method makes the assumption that:
+   *  - spacedim == meshdim (2 here).
+   *  - the point is within the quad
+   *  Quadratic elements are not supported yet.
+   *
+   *  A quadrangle can be described as 3 vectors, one point being taken as the origin.
+   *  Denoting A, B, C the three other points, any point P within the quad is written as
+   *    P = xA+ yC + xy(B-A-C)
+   *  This method solve those 2 equations (one per component) for x and y.
+   *
+
+          A------B
+          |      |
+          |      |
+          0------C
+   */
+  inline void quad_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+  {
+    double prec = 1.0e-14;
+    enum { _XX=0, _YY, _ZZ };
+
+    if(n.size() != 4)
+      throw INTERP_KERNEL::Exception("INTERP_KERNEL::quad_mapped_coords : unrecognized geometric type! Only QUAD4 supported.");
+
+    double A[2] = {n[1][_XX] - n[0][_XX],  n[1][_YY] - n[0][_YY]};
+    double B[2] = {n[2][_XX] - n[0][_XX],  n[2][_YY] - n[0][_YY]};
+    double C[2] = {n[3][_XX] - n[0][_XX],  n[3][_YY] - n[0][_YY]};
+    double N[2] = {B[_XX] - A[_XX] - C[_XX], B[_YY] - A[_YY] - C[_YY]};
+    double P[2] = {p[_XX] - n[0][_XX], p[_YY] - n[0][_YY]};
+
+    // degenerated case: a rectangle:
+    if (fabs(N[0]) < prec && fabs(N[1]) < prec)
+      {
+        double det = C[0]*A[1] -C[1]*A[0];
+        if (fabs(det) < prec)
+          throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords() has a degenerated 2x2 system!");
+        bc[0] = (P[0]*A[1]-P[1]*A[0])/det;
+        bc[1] = (P[1]*C[0]-P[0]*C[1])/det;
+        return;
+      }
+    double b,c ,a = A[1]*N[0]-A[0]*N[1];
+    bool cas1;
+    if (fabs(a) > 1.0e-14)
+      {
+        b = A[1]*C[0]+N[1]*P[0]-N[0]*P[1]-A[0]*C[1];
+        c = P[0]*C[1] - P[1]*C[0];
+        cas1 = true;
+      }
+    else
+      {
+        a = -C[1]*N[0]+C[0]*N[1];
+        b = A[1]*C[0]-N[1]*P[0]+N[0]*P[1]-A[0]*C[1];
+        c = -P[0]*A[1] + P[1]*A[0];
+        cas1 = false;
+      }
+    double delta = b*b - 4.0*a*c;
+    if (delta < 0.0)
+      throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): imaginary solutions!");
+    bc[1] = 0.5*(-b+sqrt(delta))/a;
+    if (bc[1] < -prec || bc[1] > (1.0+prec))
+      bc[1] = 0.5*(-b-sqrt(delta))/a;
+    if (bc[1] < -prec || bc[1] > (1.0+prec))
+      throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+    if (cas1)
+      {
+        double denom = C[0]+bc[1]*N[0];
+        if (fabs(denom) < prec)
+          throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+        bc[0] = (P[0]-bc[1]*A[0])/denom;
+        if (bc[0] < -prec || bc[0] > (1.0+prec))
+          throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
+      }
+    else
+      {
+        bc[0] = bc[1];
+        double denom = A[1]+bc[0]*N[1];
+        if (fabs(denom) < prec)
+          throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+        bc[1] = (P[1]-bc[0]*C[1])/denom;
+        if (bc[1] < -prec || bc[1] > (1.0+prec))
+          throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
+      }
+  }
+
+  /*!
+   * Doing as in quad_mapped_coords() would lead to a 4th order equation ... So go simpler here:
+   * orthogonal distance to each pair of parallel faces is computed. The ratio gives a number in [0,1]
+   *
+   * Conventions:
+   *   - for HEXA8, point F (5) is taken to be the origin (see med file ref connec):
+   *          0 ------ 3
+             /|       /|
+            / |      / |
+           1 ------ 2  |
+           |  |     |  |
+           |  |     |  |
+           |  4-----|- 7
+           | /      | /
+           5 ------ 6
+
+   *
+   */
+
+  inline void cuboid_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
+  {
+    double prec = 1.0e-14;
+    enum { _XX=0, _YY };
+    if (n.size() != 8)
+      throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: unrecognized geometric type! Only HEXA8 supported.");
+
+    double dx1, dx2, dy1, dy2, dz1, dz2;
+    dx1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[5],n[1]);
+    dx2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[7],n[3],n[2]);
+
+    dy1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[6],n[2]);
+    dy2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[0],n[3]);
+
+    dz1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[4],n[7]);
+    dz2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[1],n[2],n[3]);
+
+    if (dx1 < -prec || dx2 < -prec || dy1 < -prec || dy2 < -prec || dz1 < -prec || dz2 < -prec)
+      throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: point outside HEXA8");
+
+    bc[0] = dx1+dx2 < prec ? 0.5 : dx1/(dx1+dx2);
+    bc[1] = dy1+dy2 < prec ? 0.5 : dy1/(dy1+dy2);
+    bc[2] = dz1+dz2 < prec ? 0.5 : dz1/(dz1+dz2);
+  }
+
   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
   /*         calcul la surface d'un polygone.                 */
   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */

diff --git a/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.hxx b/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.hxx
new file mode 100644 (file)
index 0000000..ae70737
--- /dev/null
+++ b/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.hxx
@@ -0,0 +1,44 @@
+// Copyright (C) 2007-2016  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, or (at your option) any later version.
+//
+// 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
+//
+// Author : Adrien Bruneton (CEA/DEN)
+
+#ifndef __MappedBarycenter2DIntersectorP1P1_HXX__
+#define __MappedBarycenter2DIntersectorP1P1_HXX__
+
+#include "PlanarIntersector.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  class MappedBarycentric2DIntersectorP1P1 : public PlanarIntersector<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:
+    MappedBarycentric2DIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS, double dimCaracteristic, double md3DSurf, double minDot3DSurf, double medianPlane, double precision, int orientation);
+    void intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res);
+    int getNumberOfRowsOfResMatrix() const;
+    int getNumberOfColsOfResMatrix() const;
+  };
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.txx b/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.txx
new file mode 100644 (file)
index 0000000..65fee50
--- /dev/null
+++ b/src/INTERP_KERNEL/MappedBarycentric2DIntersectorP1P1.txx
@@ -0,0 +1,111 @@
+// Copyright (C) 2007-2016  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, or (at your option) any later version.
+//
+// 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
+//
+// Author : Adrien Bruneton (CEA/DEN)
+
+#ifndef __MappedBarycentric2DIntersectorP1P1_TXX__
+#define __MappedBarycentric2DIntersectorP1P1_TXX__
+
+#include "MappedBarycentric2DIntersectorP1P1.hxx"
+#include "PlanarIntersector.txx"
+#include "CellModel.hxx"
+
+#include "PointLocatorAlgos.txx"
+#include "MeshUtils.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  MappedBarycentric2DIntersectorP1P1<MyMeshType,MyMatrix>::MappedBarycentric2DIntersectorP1P1(const MyMeshType& meshT, const MyMeshType& meshS,
+                                                                        double dimCaracteristic, double md3DSurf, double minDot3DSurf,
+                                                                        double medianPlane, double precision, int orientation):
+    PlanarIntersector<MyMeshType,MyMatrix>(meshT,meshS,dimCaracteristic,precision,md3DSurf,minDot3DSurf,medianPlane,true,orientation,0)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  void MappedBarycentric2DIntersectorP1P1<MyMeshType,MyMatrix>::intersectCells(ConnType icellT, const std::vector<ConnType>& icellsS, MyMatrix& res)
+  {
+    std::vector<double> CoordsT;
+    PlanarIntersector<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(icellT),CoordsT);
+    int nbOfNodesT=CoordsT.size()/SPACEDIM;
+    for(typename std::vector<ConnType>::const_iterator iter=icellsS.begin();iter!=icellsS.end();iter++)
+      {
+        NormalizedCellType tS=PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getTypeOfElement(OTT<ConnType,numPol>::indFC(*iter));
+        if(tS!=NORM_QUAD4)
+          throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==2. Only QUAD4 supported !");
+        std::vector<double> CoordsS;
+        PlanarIntersector<MyMeshType,MyMatrix>::getRealSourceCoordinates(OTT<ConnType,numPol>::indFC(*iter),CoordsS);
+        std::vector<double> CoordsTTmp(CoordsT);
+        if(SPACEDIM==3)
+          PlanarIntersector<MyMeshType,MyMatrix>::projectionThis(&CoordsS[0],&CoordsTTmp[0],CoordsS.size()/SPACEDIM,nbOfNodesT);
+        const ConnType *startOfCellNodeConnT=PlanarIntersector<MyMeshType,MyMatrix>::_connectT+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexT[icellT]);
+        for(int nodeIdT=0;nodeIdT<nbOfNodesT;nodeIdT++)
+          {
+            typename MyMatrix::value_type& resRow=res[OTT<ConnType,numPol>::ind2C(startOfCellNodeConnT[nodeIdT])];
+            if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg2D(&CoordsTTmp[nodeIdT*SPACEDIM],&CoordsS[0],4,PlanarIntersector<MyMeshType,MyMatrix>::_precision) )
+              {
+                double mco[2];  // mapped coordinates in the quad4
+                std::vector<const double*> coo(4);
+                coo[0]=&CoordsS[0]; coo[1]=&CoordsS[SPACEDIM]; coo[2]=&CoordsS[2*SPACEDIM]; coo[3]=&CoordsS[3*SPACEDIM];
+                quad_mapped_coords(coo,&CoordsTTmp[nodeIdT*SPACEDIM],mco);
+
+                // Now use the form function of the QUAD4 to map the field values
+                double resLoc[4];
+                // See QUAD4 standard connectivity and cuboid_mapped_coords() convention:
+                resLoc[0] = (1.-mco[0]) * (1.-mco[1]);
+                resLoc[1] = (1.-mco[0]) *   mco[1]   ;
+                resLoc[2] =  mco[0]     *   mco[1]   ;
+                resLoc[3] =  mco[0]     * (1.-mco[1]);
+
+                const ConnType *startOfCellNodeConnS=PlanarIntersector<MyMeshType,MyMatrix>::_connectS+OTT<ConnType,numPol>::conn2C(PlanarIntersector<MyMeshType,MyMatrix>::_connIndexS[*iter]);
+                for(int nodeIdS=0;nodeIdS<4;nodeIdS++)
+                  {
+                    if(fabs(resLoc[nodeIdS])>PlanarIntersector<MyMeshType,MyMatrix>::_precision)
+                      {
+                        ConnType curNodeSInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConnS[nodeIdS]);
+                        typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
+                        if(iterRes==resRow.end())
+                          resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(curNodeSInCmode),resLoc[nodeIdS]));
+                        else
+                          {
+                            double val=(*iterRes).second+resLoc[nodeIdS];
+                            resRow.erase(OTT<ConnType,numPol>::indFC(curNodeSInCmode));
+                            resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(curNodeSInCmode),val));
+                          }
+                      }
+                  }
+              }
+          }
+      }
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int MappedBarycentric2DIntersectorP1P1<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
+  {
+    return PlanarIntersector<MyMeshType,MyMatrix>::_meshT.getNumberOfNodes();
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int MappedBarycentric2DIntersectorP1P1<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
+  {
+    return PlanarIntersector<MyMeshType,MyMatrix>::_meshS.getNumberOfNodes();
+  }
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.hxx b/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.hxx
new file mode 100644 (file)
index 0000000..0085d11
--- /dev/null
+++ b/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.hxx
@@ -0,0 +1,46 @@
+// Copyright (C) 2007-2016  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, or (at your option) any later version.
+//
+// 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
+//
+// Author : Adrien Bruneton (CEA/DEN)
+
+#ifndef __MappedBarycentric3DIntersectorP1P1_HXX__
+#define __MappedBarycentric3DIntersectorP1P1_HXX__
+
+#include "Intersector3DP1P1.hxx"
+#include "NormalizedUnstructuredMesh.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  class MappedBarycentric3DIntersectorP1P1 : public Intersector3DP1P1<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:
+    MappedBarycentric3DIntersectorP1P1(const MyMeshType& targetMesh, const MyMeshType& srcMesh, double precision);
+    ~MappedBarycentric3DIntersectorP1P1();
+    void intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res);
+  protected:
+    double _precision;
+  };
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.txx b/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.txx
new file mode 100644 (file)
index 0000000..018d5b7
--- /dev/null
+++ b/src/INTERP_KERNEL/MappedBarycentric3DIntersectorP1P1.txx
@@ -0,0 +1,113 @@
+// Copyright (C) 2007-2016  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, or (at your option) any later version.
+//
+// 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
+//
+// Author : Adrien Bruneton (CEA/DEN)
+
+#ifndef __MAPPEDBARYCENTRIC3DINTERSECTORP1P1_TXX__
+#define __MAPPEDBARYCENTRIC3DINTERSECTORP1P1_TXX__
+
+#include "MappedBarycentric3DIntersectorP1P1.hxx"
+#include "Intersector3DP1P1.txx"
+#include "MeshUtils.hxx"
+
+namespace INTERP_KERNEL
+{
+
+  /**
+   * Constructor creating object from target cell global number 
+   * 
+   * @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>
+  MappedBarycentric3DIntersectorP1P1<MyMeshType,MyMatrix>::MappedBarycentric3DIntersectorP1P1(const MyMeshType& targetMesh, const MyMeshType& srcMesh, double precision):
+    Intersector3DP1P1<MyMeshType,MyMatrix>(targetMesh,srcMesh),_precision(precision)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  MappedBarycentric3DIntersectorP1P1<MyMeshType,MyMatrix>::~MappedBarycentric3DIntersectorP1P1()
+  {
+  }
+
+  /**
+   * @param targetCell in C mode.
+   * @param srcCells in C mode.
+   */
+  template<class MyMeshType, class MyMatrix>
+  void MappedBarycentric3DIntersectorP1P1<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
+  {
+    std::vector<double> CoordsT;
+    const ConnType *startOfCellNodeConnT=Intersector3DP1P1<MyMeshType,MyMatrix>::getStartConnOfTargetCell(targetCell);
+    Intersector3DP1P1<MyMeshType,MyMatrix>::getRealTargetCoordinates(OTT<ConnType,numPol>::indFC(targetCell),CoordsT);
+    int nbOfNodesT=CoordsT.size()/SPACEDIM;
+    const double *coordsS=Intersector3DP1P1<MyMeshType,MyMatrix>::_src_mesh.getCoordinatesPtr();
+    for(int nodeIdT=0;nodeIdT<nbOfNodesT;nodeIdT++)
+      {
+        typename MyMatrix::value_type& resRow=res[OTT<ConnType,numPol>::ind2C(startOfCellNodeConnT[nodeIdT])];
+        if(!resRow.empty())
+          continue;
+        for(typename std::vector<ConnType>::const_iterator iterCellS=srcCells.begin();iterCellS!=srcCells.end();iterCellS++)
+          {
+            NormalizedCellType tS=Intersector3DP1P1<MyMeshType,MyMatrix>::_src_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(*iterCellS));
+            if(tS!=NORM_HEXA8)
+              throw INTERP_KERNEL::Exception("Invalid source cell detected for meshdim==3. Only HEXA8 supported !");
+            const CellModel& cmTypeS=CellModel::GetCellModel(tS);
+            //
+            std::vector<ConnType> connOfCurCellS;
+            Intersector3DP1P1<MyMeshType,MyMatrix>::getConnOfSourceCell(OTT<ConnType,numPol>::indFC(*iterCellS),connOfCurCellS);
+            if( PointLocatorAlgos<MyMeshType>::isElementContainsPointAlg3D(&CoordsT[nodeIdT*SPACEDIM],&connOfCurCellS[0],connOfCurCellS.size(),coordsS,cmTypeS,_precision) )
+              {
+                double mco[3];  // mapped coordinates in the hexa8
+                std::vector<double> localCoordsS;
+                Intersector3DP1P1<MyMeshType,MyMatrix>::getRealSourceCoordinates(OTT<ConnType,numPol>::indFC(*iterCellS),localCoordsS);
+                std::vector<const double*> coo(8);
+                coo[0]=&localCoordsS[0]; coo[1]=&localCoordsS[3]; coo[2]=&localCoordsS[6]; coo[3]=&localCoordsS[9];
+                coo[4]=&localCoordsS[12]; coo[5]=&localCoordsS[15]; coo[6]=&localCoordsS[18]; coo[7]=&localCoordsS[21];
+                cuboid_mapped_coords(coo,&CoordsT[nodeIdT*SPACEDIM],mco);
+
+                // Now use the form function of the HEXA8 to map the field values
+                double resLoc[8];
+                // See HEXA8 standard connectivity and cuboid_mapped_coords() convention:
+                resLoc[5] = (1.-mco[0]) * (1.-mco[1]) * (1.-mco[2]);
+                resLoc[6] =  mco[0]     * (1.-mco[1]) * (1.-mco[2]);
+                resLoc[7] =  mco[0]     *   mco[1]    * (1.-mco[2]);
+                resLoc[4] = (1.-mco[0]) *   mco[1]    * (1.-mco[2]);
+
+                resLoc[1] = (1.-mco[0]) * (1.-mco[1]) * mco[2];
+                resLoc[2] =  mco[0]     * (1.-mco[1]) * mco[2];
+                resLoc[3] =  mco[0]     *   mco[1]    * mco[2];
+                resLoc[0] = (1.-mco[0]) *   mco[1]    * mco[2];
+
+                const ConnType *startOfCellNodeConnS=Intersector3DP1P1<MyMeshType,MyMatrix>::getStartConnOfSourceCell(*iterCellS);
+                for(int nodeIdS=0;nodeIdS<8;nodeIdS++)
+                  {
+                    if(fabs(resLoc[nodeIdS])>_precision)
+                      {
+                        ConnType curNodeSInCmode=OTT<ConnType,numPol>::coo2C(startOfCellNodeConnS[nodeIdS]);
+                        resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(curNodeSInCmode),resLoc[nodeIdS]));
+                      }
+                  }
+              }
+          }
+      }
+  }
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/VolSurfUser.cxx b/src/INTERP_KERNEL/VolSurfUser.cxx
index 3aa6c9ae0db5d4ea3227eb190180cac8bf31502d..2d06cdde06aa8e3720b14925de3421e44c065fb9 100644 (file)
--- a/src/INTERP_KERNEL/VolSurfUser.cxx
+++ b/src/INTERP_KERNEL/VolSurfUser.cxx
@@ -28,6 +28,25 @@
 
 namespace INTERP_KERNEL
 {
+  /* Orthogonal distance from a point to a plane defined by three points p1, p2, p3.
+   * Returns a signed distance, the normal of the plane being defined by (p1-p2)x(p3-p2)
+   */
+  double OrthoDistanceFromPtToPlaneInSpaceDim3(const double *p, const double *p1, const double *p2, const double *p3)
+  {
+    double prec = 1.0e-14;
+    double T[2][3] = {{p1[0] - p2[0], p1[1] - p2[1], p1[2] - p2[2]},
+                      {p3[0] - p2[0], p3[1] - p2[1], p3[2] - p2[2]}};
+    double N[3] = {T[0][1]*T[1][2]-T[0][2]*T[1][1],
+                   T[0][2]*T[1][0]-T[0][0]*T[1][2],
+                   T[0][0]*T[1][1]-T[0][1]*T[1][0]};
+
+    double norm2 = N[0]*N[0] + N[1]*N[1] + N[2]*N[2];
+    if (norm2 < prec)
+      throw INTERP_KERNEL::Exception("OrthoDistanceFromPtToPlaneInSpaceDim3: degenerated normal vector!");
+    double num = N[0]*(p[0]-p1[0]) + N[1]*(p[1]-p1[1]) + N[2]*(p[2]-p1[2]);
+    return num/sqrt(norm2);
+  }
+
   double SquareDistanceFromPtToSegInSpaceDim2(const double *pt, const double *pt0Seg2, const double *pt1Seg2, std::size_t &nbOfHint)
   {
     double dx=pt1Seg2[0]-pt0Seg2[0],dy=pt1Seg2[1]-pt0Seg2[1];

diff --git a/src/INTERP_KERNEL/VolSurfUser.hxx b/src/INTERP_KERNEL/VolSurfUser.hxx
index ece473fcf07a30e1b46e0eda60a2b44d2af37107..25938345bfae390203dc4879a828e21b09a83514 100644 (file)
--- a/src/INTERP_KERNEL/VolSurfUser.hxx
+++ b/src/INTERP_KERNEL/VolSurfUser.hxx
@@ -39,6 +39,8 @@ namespace INTERP_KERNEL
   template<class ConnType, NumberingPolicy numPolConn>
   void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res);
 
+  double INTERPKERNEL_EXPORT OrthoDistanceFromPtToPlaneInSpaceDim3(const double *p, const double *p1, const double *p2, const double *p3);
+
   double INTERPKERNEL_EXPORT SquareDistanceFromPtToSegInSpaceDim2(const double *pt, const double *pt0Seg2, const double *pt1Seg2, std::size_t &nbOfHint);
 
   double INTERPKERNEL_EXPORT DistanceFromPtToTriInSpaceDim3(const double *pt, const double *pt0Tri3, const double *pt1Tri3, const double *pt2Tri3);

diff --git a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx
index add2bfea276589739559fb6d0db872b95f502281..851685eaae4fc49f86dafa4fd1243aba79f952a6 100644 (file)
--- a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx
+++ b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.cxx
@@ -342,5 +342,148 @@ namespace INTERP_TEST
     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[0], p2[0], 1e-12);
     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[1], p2[1], 1e-12);
     CPPUNIT_ASSERT_DOUBLES_EQUAL( p[2], p2[2], 1e-12);
-  }  
+  }
+
+  /* Conventions:
+  *   - for HEXA8, point 5 is taken to be the origin (see med file ref connec):
+  *          0 ------ 3
+            /|       /|
+           / |      / |
+          1 ------ 2  |
+          |  |     |  |
+          |  |     |  |
+          |  4-----|- 7
+          | /      | /
+          5 ------ 6
+   */
+  void UnitTetraIntersectionBaryTest::test_cuboid_mapped_coords_3D()
+  {
+    double nodes[8][3] = { { 0.0, 2.0, 4.0 }, //0
+                           { 0.0, 0.0, 4.0 },
+                           { 1.0, 0.0, 4.0 },
+                           { 1.0, 2.0, 4.0 },
+                           { 0.0, 2.0, 0.0 }, // 4
+                           { 0.0, 0.0, 0.0 },
+                           { 1.0, 0.0, 0.0 },
+                           { 1.0, 2.0, 0.0 }
+    };
+    // Translate cube:
+    for (int i=0; i < 8; ++i)
+      for (int j=0; j < 3; ++j)
+        nodes[i][j] += 15.0;
+
+    std::vector<const double*> n (8);
+    for (int i=0; i<8; i++)
+      n[i] = &nodes[i][0];
+
+    {
+        // middle point
+        double p[3] = { 15.5, 16.0, 17.0 }, bc[3];
+        cuboid_mapped_coords(n, p, bc);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
+        CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[2], 1e-12);
+    }
+    {
+      // point 1
+      double p[3] = { 15.0, 15.0, 19.0 }, bc[3];
+      cuboid_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
+    }
+    {
+      // point 7
+      double p[3] = { 16.0, 17.0, 15.0 }, bc[3];
+      cuboid_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[2], 1e-12);
+    }
+    {
+      // point 3
+      double p[3] = { 16.0, 17.0, 19.0 }, bc[3];
+      cuboid_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[2], 1e-12);
+    }
+    {
+      // point outside
+      double p[3] = { 2.0, 16.0, 18.0 }, bc[3];
+      CPPUNIT_ASSERT_THROW(cuboid_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
+    }
+
+  }
+
+  /* Convention
+      - for QUAD4, point 0 is taken to be the origin (again see med file ref connec):
+
+         1------2
+         |      |
+         |      |
+         0------3
+  */
+  void UnitTetraIntersectionBaryTest::test_quad_mapped_coords_2D()
+  {
+
+    double nodes[4][2] = { { 0.0, 0.0 },
+                           { 0.0, 1.0 },
+                           { 2.0, 3.0 },
+                           { 1.0, 0.0 } };
+
+    // Translate quad4:
+    for (int i=0; i < 4; ++i)
+      for (int j=0; j < 2; ++j)
+        nodes[i][j] += 15.0;
+
+    std::vector<const double*> n (4);
+    for (int i=0; i<4; i++)
+      n[i] = &nodes[i][0];
+
+    {
+      // middle point
+      double p[2] = { 15.75, 16.0 }, bc[2];
+      quad_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[1], 1e-12);
+    }
+
+    {
+      // middle point of seg
+      double p[2] = { 15.5, 15.0 }, bc[2];
+      quad_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
+    }
+
+    {
+      // point 1
+      double p[2] = { 15.0, 16.0 }, bc[2];
+      quad_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
+    }
+    {
+      // point 2
+      double p[2] = { 17.0, 18.0 }, bc[2];
+      quad_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[1], 1e-12);
+    }
+    {
+      // point 3
+      double p[2] = { 16.0, 15.0 }, bc[2];
+      quad_mapped_coords(n, p, bc);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, bc[0], 1e-12);
+      CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, bc[1], 1e-12);
+    }
+    {
+      // point outside
+      double p[2] = { 18.0, 18.0 }, bc[2];
+      CPPUNIT_ASSERT_THROW(quad_mapped_coords(n, p, bc), INTERP_KERNEL::Exception);
+    }
+  }
+
+
 }

diff --git a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx
index 952481e5f68e3979759a28ce268430b3f0d45a73..3657620b1902f3b30198afd956568156a4e2c8ca 100644 (file)
--- a/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx
+++ b/src/INTERP_KERNELTest/UnitTetraIntersectionBaryTest.hxx
@@ -52,6 +52,8 @@ namespace INTERP_TEST
     CPPUNIT_TEST( test_UnitTetraIntersectionBary_11 );
     CPPUNIT_TEST( test_TetraAffineTransform_reverseApply );
     CPPUNIT_TEST( test_barycentric_coords );
+    CPPUNIT_TEST( test_cuboid_mapped_coords_3D );
+    CPPUNIT_TEST( test_quad_mapped_coords_2D );
     CPPUNIT_TEST_SUITE_END();
   public:
     void test_UnitTetraIntersectionBary_1();
@@ -69,6 +71,8 @@ namespace INTERP_TEST
     void test_UnitTetraIntersectionBary_13();
     void test_TetraAffineTransform_reverseApply();
     void test_barycentric_coords();
+    void test_cuboid_mapped_coords_3D();
+    void test_quad_mapped_coords_2D();
   };
 }
 

diff --git a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
index 1398ffebe88a993e23b4294c3acc0bdc845beec7..5b4b2261fba941c1205da584e5857d870052af24 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingRemapperTest.py
@@ -555,7 +555,7 @@ class MEDCouplingBasicsTest(unittest.TestCase):
     def testSwig2BarycentricP1P13D_1(self):
         sCoo=DataArrayDouble([0.313,0.00218,6.90489,0.313,0.10692667,6.90489,0.313,0.10692667,6.96790167,0.313,0.00218,6.9773125,0.313,0.21167333,6.90489,0.313,0.21167333,6.95849083,0.313,0.31642,6.90489,0.313,0.31642,6.94908,0.313,0.09383333,7.04891667,0.313,0.00218,7.049735,0.313,0.18548667,7.04809833,0.313,0.27714,7.04728,0.313,0.05782667,7.133205,0.313,0.00218,7.1221575,0.313,0.11347333,7.1442525,0.313,0.16912,7.1553,0.313,0.02509333,7.19458,0.313,0.00218,7.19458,0.313,0.04800667,7.19458,0.313,0.07092,7.19458,0.31005609,0.00218,6.90460005,0.31005609,0.10692667,6.90460005,0.29776312,0.10692667,6.96640097,0.29592716,0.00218,6.97563097,0.31005609,0.21167333,6.90460005,0.29959908,0.21167333,6.95717096,0.31005609,0.31642,6.90460005,0.30143505,0.31642,6.94794095,0.28195788,0.09383333,7.04585928,0.28179823,0.00218,7.04666189,0.28211753,0.18548667,7.04505668,0.28227718,0.27714,7.04425407,0.26551404,0.05782667,7.12852804,0.2676693,0.00218,7.11769282,0.26335878,0.11347333,7.13936327,0.26120352,0.16912,7.15019849,0.25354037,0.02509333,7.18872374,0.25354037,0.00218,7.18872374,0.25354037,0.04800667,7.18872374,0.25354037,0.07092,7.18872374,0.30722531,0.00218,6.90374134,0.30722531,0.10692667,6.90374134,0.28311179,0.10692667,6.96195653,0.27951042,0.00218,6.97065101,0.30722531,0.21167333,6.90374134,0.28671316,0.21167333,6.95326205,0.30722531,0.31642,6.90374134,0.29031453,0.31642,6.94456758,0.25210869,0.09383333,7.03680463,0.25179553,0.00218,7.03756067,0.25242185,0.18548667,7.03604859,0.25273501,0.27714,7.03529255,0.21985294,0.05782667,7.1146769,0.22408063,0.00218,7.10447034,0.21562524,0.11347333,7.12488346,0.21139755,0.16912,7.13509002,0.19636574,0.02509333,7.17138,0.19636574,0.00218,7.17138,0.19636574,0.04800667,7.17138,0.19636574,0.07092,7.17138,0.30461645,0.00218,6.90234688,0.30461645,0.10692667,6.90234688,0.26960904,0.10692667,6.95473916,0.26438066,0.00218,6.96256398,0.30461645,0.21167333,6.90234688,0.27483742,0.21167333,6.94691434,0.30461645,0.31642,6.90234688,0.2800658,0.31642,6.93908952,0.22459952,0.09383333,7.02210067,0.22414487,0.00218,7.02278109,0.22505416,0.18548667,7.02142025,0.2255088,0.27714,7.02073983,0.17777143,0.05782667,7.09218386,0.18390909,0.00218,7.0829982,0.17163377,0.11347333,7.10136952,0.1654961,0.16912,7.11055518,0.1436733,0.02509333,7.14321531,0.1436733,0.00218,7.14321531,0.1436733,0.04800667,7.14321531,0.1436733,0.07092,7.14321531,0.30232976,0.00218,6.90047024,0.30232976,0.10692667,6.90047024,0.25777378,0.10692667,6.94502622,0.25111932,0.00218,6.95168068,0.30232976,0.21167333,6.90047024,0.26442825,0.21167333,6.93837175,0.30232976,0.31642,6.90047024,0.27108271,0.31642,6.93171729,0.20048753,0.09383333,7.00231247,0.19990888,0.00218,7.00289112,0.20106618,0.18548667,7.00173382,0.20164482,0.27714,7.00115518,0.14088667,0.05782667,7.06191333,0.14869844,0.00218,7.05410156,0.13307491,0.11347333,7.06972509,0.12526315,0.16912,7.07753685,0.097488,0.02509333,7.105312,0.097488,0.00218,7.105312,0.097488,0.04800667,7.105312,0.097488,0.07092,7.105312,0.30045312,0.00218,6.89818355,0.30045312,0.10692667,6.89818355,0.24806084,0.10692667,6.93319096,0.24023602,0.00218,6.93841934,0.30045312,0.21167333,6.89818355,0.25588566,0.21167333,6.92796258,0.30045312,0.31642,6.89818355,0.26371048,0.31642,6.9227342,0.18069933,0.09383333,6.97820048,0.18001891,0.00218,6.97865513,0.18137975,0.18548667,6.97774584,0.18206017,0.27714,6.9772912,0.11061614,0.05782667,7.02502857,0.1198018,0.00218,7.01889091,0.10143048,0.11347333,7.03116623,0.09224482,0.16912,7.0373039,0.05958469,0.02509333,7.0591267,0.05958469,0.00218,7.0591267,0.05958469,0.04800667,7.0591267,0.05958469,0.07092,7.0591267,0.29905866,0.00218,6.89557469,0.29905866,0.10692667,6.89557469,0.24084347,0.10692667,6.91968821,0.23214899,0.00218,6.92328958,0.29905866,0.21167333,6.89557469,0.24953795,0.21167333,6.91608684,0.29905866,0.31642,6.89557469,0.25823242,0.31642,6.91248547,0.16599537,0.09383333,6.95069131,0.16523933,0.00218,6.95100447,0.16675141,0.18548667,6.95037815,0.16750745,0.27714,6.95006499,0.0881231,0.05782667,6.98294706,0.09832966,0.00218,6.97871937,0.07791654,0.11347333,6.98717476,0.06770998,0.16912,6.99140245,0.03142,0.02509333,7.00643426,0.03142,0.00218,7.00643426,0.03142,0.04800667,7.00643426,0.03142,0.07092,7.00643426,0.29819995,0.00218,6.89274391,0.29819995,0.10692667,6.89274391,0.23639903,0.10692667,6.90503688,0.22716903,0.00218,6.90687284,0.29819995,0.21167333,6.89274391,0.24562904,0.21167333,6.90320092,0.29819995,0.31642,6.89274391,0.25485905,0.31642,6.90136495,0.15694072,0.09383333,6.92084212,0.15613811,0.00218,6.92100177,0.15774332,0.18548667,6.92068247,0.15854593,0.27714,6.92052282,0.07427196,0.05782667,6.93728596,0.08510718,0.00218,6.9351307,0.06343673,0.11347333,6.93944122,0.05260151,0.16912,6.94159648,0.01407626,0.02509333,6.94925963,0.01407626,0.00218,6.94925963,0.01407626,0.04800667,6.94925963,0.01407626,0.07092,6.94925963,0.29792818,0.00218,6.89054043,0.29792818,0.10692667,6.89054043,0.23499241,0.10692667,6.89363227,0.22559291,0.00218,6.89409403,0.29792818,0.21167333,6.89054043,0.24439191,0.21167333,6.8931705,0.29792818,0.31642,6.89054043,0.25379141,0.31642,6.89270873,0.154075,0.09383333,6.89760748,0.15325765,0.00218,6.89764764,0.15489234,0.18548667,6.89756733,0.15570969,0.27714,6.89752718,0.06988819,0.05782667,6.90174332,0.08092238,0.00218,6.90120124,0.058854,0.11347333,6.90228539,0.04781981,0.16912,6.90282747,0.00858712,0.02509333,6.90475485,0.00858712,0.00218,6.90475485,0.00858712,0.04800667,6.90475485,0.00858712,0.07092,6.90475485,0.29791,0.00218,6.820902,0.29791,0.10692667,6.820902,0.23489833,0.10692667,6.820902,0.2254875,0.00218,6.820902,0.29791,0.21167333,6.820902,0.24430917,0.21167333,6.820902,0.29791,0.31642,6.820902,0.25372,0.31642,6.820902,0.15388333,0.09383333,6.820902,0.153065,0.00218,6.820902,0.15470167,0.18548667,6.820902,0.15552,0.27714,6.820902,0.069595,0.05782667,6.820902,0.0806425,0.00218,6.820902,0.0585475,0.11347333,6.820902,0.0475,0.16912,6.820902,0.00822,0.02509333,6.820902,0.00822,0.00218,6.820902,0.00822,0.04800667,6.820902,0.00822,0.07092,6.820902],200,3)
         sConn=DataArrayInt([0,1,2,3,20,21,22,23,1,4,5,2,21,24,25,22,4,6,7,5,24,26,27,25,3,2,8,9,23,22,28,29,2,5,10,8,22,25,30,28,5,7,11,10,25,27,31,30,9,8,12,13,29,28,32,33,8,10,14,12,28,30,34,32,10,11,15,14,30,31,35,34,13,12,16,17,33,32,36,37,12,14,18,16,32,34,38,36,14,15,19,18,34,35,39,38,20,21,22,23,40,41,42,43,21,24,25,22,41,44,45,42,24,26,27,25,44,46,47,45,23,22,28,29,43,42,48,49,22,25,30,28,42,45,50,48,25,27,31,30,45,47,51,50,29,28,32,33,49,48,52,53,28,30,34,32,48,50,54,52,30,31,35,34,50,51,55,54,33,32,36,37,53,52,56,57,32,34,38,36,52,54,58,56,34,35,39,38,54,55,59,58,40,41,42,43,60,61,62,63,41,44,45,42,61,64,65,62,44,46,47,45,64,66,67,65,43,42,48,49,63,62,68,69,42,45,50,48,62,65,70,68,45,47,51,50,65,67,71,70,49,48,52,53,69,68,72,73,48,50,54,52,68,70,74,72,50,51,55,54,70,71,75,74,53,52,56,57,73,72,76,77,52,54,58,56,72,74,78,76,54,55,59,58,74,75,79,78,60,61,62,63,80,81,82,83,61,64,65,62,81,84,85,82,64,66,67,65,84,86,87,85,63,62,68,69,83,82,88,89,62,65,70,68,82,85,90,88,65,67,71,70,85,87,91,90,69,68,72,73,89,88,92,93,68,70,74,72,88,90,94,92,70,71,75,74,90,91,95,94,73,72,76,77,93,92,96,97,72,74,78,76,92,94,98,96,74,75,79,78,94,95,99,98,80,81,82,83,100,101,102,103,81,84,85,82,101,104,105,102,84,86,87,85,104,106,107,105,83,82,88,89,103,102,108,109,82,85,90,88,102,105,110,108,85,87,91,90,105,107,111,110,89,88,92,93,109,108,112,113,88,90,94,92,108,110,114,112,90,91,95,94,110,111,115,114,93,92,96,97,113,112,116,117,92,94,98,96,112,114,118,116,94,95,99,98,114,115,119,118,100,101,102,103,120,121,122,123,101,104,105,102,121,124,125,122,104,106,107,105,124,126,127,125,103,102,108,109,123,122,128,129,102,105,110,108,122,125,130,128,105,107,111,110,125,127,131,130,109,108,112,113,129,128,132,133,108,110,114,112,128,130,134,132,110,111,115,114,130,131,135,134,113,112,116,117,133,132,136,137,112,114,118,116,132,134,138,136,114,115,119,118,134,135,139,138,120,121,122,123,140,141,142,143,121,124,125,122,141,144,145,142,124,126,127,125,144,146,147,145,123,122,128,129,143,142,148,149,122,125,130,128,142,145,150,148,125,127,131,130,145,147,151,150,129,128,132,133,149,148,152,153,128,130,134,132,148,150,154,152,130,131,135,134,150,151,155,154,133,132,136,137,153,152,156,157,132,134,138,136,152,154,158,156,134,135,139,138,154,155,159,158,140,141,142,143,160,161,162,163,141,144,145,142,161,164,165,162,144,146,147,145,164,166,167,165,143,142,148,149,163,162,168,169,142,145,150,148,162,165,170,168,145,147,151,150,165,167,171,170,149,148,152,153,169,168,172,173,148,150,154,152,168,170,174,172,150,151,155,154,170,171,175,174,153,152,156,157,173,172,176,177,152,154,158,156,172,174,178,176,154,155,159,158,174,175,179,178,160,161,162,163,180,181,182,183,161,164,165,162,181,184,185,182,164,166,167,165,184,186,187,185,163,162,168,169,183,182,188,189,162,165,170,168,182,185,190,188,165,167,171,170,185,187,191,190,169,168,172,173,189,188,192,193,168,170,174,172,188,190,194,192,170,171,175,174,190,191,195,194,173,172,176,177,193,192,196,197,172,174,178,176,192,194,198,196,174,175,179,178,194,195,199,198])
-        s=MEDCoupling1SGTUMesh("target",NORM_HEXA8) ; s.setCoords(sCoo)
+        s=MEDCoupling1SGTUMesh("source",NORM_HEXA8) ; s.setCoords(sCoo)
         s.setNodalConnectivity(sConn)
         #
         tCoo=DataArrayDouble([0.328,0.012,6.8598,0.328,0.168320184237353,6.8598,0.328,0.324640368474706,6.8598,0.328,0.0,6.8598,0.298,0.012,6.8598,0.1565,0.012,6.8598,0.180205346493166,0.144794653506834,6.8598,0.298,0.168320184237353,6.8598,0.0,0.012,6.8598,0.0916755774886107,0.233324422511389,6.8598,0.298,0.324640368474706,6.8598,0.298,0.0,6.8598,0.1565,0.0,6.8598,0.0,0.0,6.8598,0.328,0.012,7.2298,0.328,0.168320184237353,7.2298,0.328,0.324640368474706,7.2298,0.328,0.0,7.2298,0.298,0.012,7.2298,0.1565,0.012,7.2298,0.180205346493166,0.144794653506834,7.2298,0.298,0.168320184237353,7.2298,0.0,0.012,7.2298,0.0916755774886107,0.233324422511389,7.2298,0.298,0.324640368474706,7.2298,0.298,0.0,7.2298,0.1565,0.0,7.2298,0.0,0.0,7.2298],28,3)
@@ -587,6 +587,113 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertAlmostEqual(0.3521445110626687 ,m[6][170],12)
         pass
 
+    def testSwig2MappedBarycentricP1P12D_1(self):
+        """ Testing mapped barycentric P1P1 projection
+        (uses analytical mapping from square to arbitrary convex quadrangle)
+        """
+        n = 5
+        sCoo = DataArrayDouble(n,1)
+        sCoo.iota(0.0);     sCoo /= float(n-1)
+        m = MEDCouplingCMesh("target")
+        m.setCoordsAt(0, sCoo)
+        m.setCoordsAt(1, sCoo)
+        tgt = m.buildUnstructured()
+        coo = tgt.getCoords()
+        orig = coo.deepCopy();   orig[:,0] = 10.0; orig[:,1] = 15.0
+        pt_a = coo.deepCopy();   pt_a[:,0] = -0.3; pt_a[:,1] = 1.0
+        pt_b = coo.deepCopy();   pt_b[:,0] = 2.0;  pt_b[:,1] = 3.0
+        pt_c = coo.deepCopy();   pt_c[:,0] = 1.0;  pt_c[:,1] = 0.0
+        # P = x*C+y*A + xy(B-A-C) + ORIGIN
+        coo2 = coo[:,0]*pt_c + coo[:, 1]*pt_a + coo[:, 0]*coo[:, 1]*(pt_b - pt_a - pt_c) + orig
+         
+        tgt.setCoords(coo2)
+      
+        sCoo = DataArrayDouble([0.0,0.0,  -0.3,1.0,  2.0,3.0,  1.0,0.0],4,2)
+        sCoo[:,0] += 10.0;  sCoo[:,1] += 15.0;   
+        sConn = DataArrayInt([0,1,2,3])
+        s = MEDCoupling1SGTUMesh("source",NORM_QUAD4) ; s.setCoords(sCoo)
+        s.setNodalConnectivity(sConn)
+        #
+        aRemapper=MEDCouplingRemapper()
+        aRemapper.setPrecision(1e-12)
+        aRemapper.setIntersectionType(MappedBarycentric)
+        self.assertEqual(aRemapper.prepare(s,tgt,'P1P1'),1)
+        srcField = MEDCouplingFieldDouble(ON_NODES, ONE_TIME)
+        srcField.setNature(IntensiveMaximum)
+        srcField.setMesh(s); srcField.setName("field")
+        srcField.setArray(DataArrayDouble([1.0,2.0,3.0,4.0]))
+        tgtF = aRemapper.transferField(srcField, 1e+300)
+        ref = [1.0, 1.75, 2.5, 3.25, 4.0, 1.25, 1.875, 2.5, 3.125, 3.75, 1.5, 2.0, 2.5, 3.0, 3.5, 1.75, 
+         2.125, 2.5, 2.875, 3.25, 2.0, 2.25, 2.5, 2.75, 3.0]
+        val = tgtF.getArray().getValues()
+        for i, ref_v in enumerate(ref):
+            self.assertAlmostEqual(ref_v, val[i])        
+        pass
+
+    def testSwig2MappedBarycentricP1P13_1(self):
+        """ Testing mapped barycentric P1P1 projection in 3D (uses orthogonal distances to 
+        HEXA8 faces).
+        Convention:
+              0 ------ 3
+             /|       /|
+            / |      / |
+           1 ------ 2  |
+           |  |     |  |
+           |  |     |  |
+           |  4-----|- 7
+           | /      | /
+           5 ------ 6
+        """
+        n = 5
+        sCoo = DataArrayDouble(n,1)
+        sCoo.iota(0.0)
+        sCoo /= float(n-1)
+        m = MEDCouplingCMesh("target")
+        m.setCoordsAt(0, sCoo)
+        m.setCoordsAt(1, sCoo)
+        m.setCoordsAt(2, sCoo)
+        tgt = m.buildUnstructured()
+        coo = tgt.getCoords()
+        pt_0 = coo.deepCopy(); pt_0[:,0] = -0.3; pt_0[:,1] = 1.0; pt_0[:,2] = 1.0
+        pt_1 = coo.deepCopy(); pt_1[:,0] = 0.0; pt_1[:,1] = 0.0; pt_1[:,2] = 1.0
+        pt_2 = coo.deepCopy(); pt_2[:,0] = 1.0; pt_2[:,1] = 0.0; pt_2[:,2] = 1.0
+        pt_3 = coo.deepCopy(); pt_3[:,0] = 2.0; pt_3[:,1] = 3.0; pt_3[:,2] = 1.0
+        
+        pt_4 = coo.deepCopy(); pt_4[:,0] = -0.3; pt_4[:,1] = 1.0; pt_4[:,2] = 0.0
+        orig = coo.deepCopy(); orig[:,0] = 10.0; orig[:,1] = 15.0; orig[:,2] = 20.0
+        pt_6 = coo.deepCopy(); pt_6[:,0] = 1.0; pt_6[:,1] = 0.0; pt_6[:,2] = 0.0 
+        pt_7 = coo.deepCopy(); pt_7[:,0] = 2.0; pt_7[:,1] = 3.0; pt_7[:,2] = 0.0
+        # P = x*p6 + y*p4 + z*p1 + xy*(p7-p6-p4) + xz*(p2-p1-p6) + yz*(p0-p4-p1) + xyz(p3-p7-p2-p0+p1+p6+p4)
+        x,y,z = coo[:,0],coo[:,1],coo[:,2]
+        coo2 = x*pt_6 + y*pt_4 + z*pt_1 + \
+               x*y*(pt_7 - pt_6 - pt_4) + x*z*(pt_2 - pt_1 - pt_6) + y*z*(pt_0 - pt_4 - pt_1) + \
+               x*y*z*(pt_3 - pt_7 - pt_2 - pt_0 + pt_6 + pt_1 + pt_4) + orig
+        tgt.setCoords(coo2)
+      
+        sCoo = DataArrayDouble([-0.3,1.0,1.0,  0.0,0.0,1.0,  1.0,0.0,1.0,  2.0,3.0,1.0, 
+                                -0.3,1.0,0.0,  0.0,0.0,0.0,  1.0,0.0,0.0,  2.0,3.0,0.0,],8,3)
+        sCoo[:, 0] += 10.0; sCoo[:, 1] += 15.0; sCoo[:, 2] += 20.0;  
+        sConn = DataArrayInt([0,1,2,3,4, 5,6,7])
+        s = MEDCoupling1SGTUMesh("source",NORM_HEXA8) ; s.setCoords(sCoo)
+        s.setNodalConnectivity(sConn)
+        #
+        aRemapper=MEDCouplingRemapper()
+        aRemapper.setPrecision(1e-12)
+        aRemapper.setIntersectionType(MappedBarycentric)
+        self.assertEqual(aRemapper.prepare(s,tgt,'P1P1'),1)
+        srcField = MEDCouplingFieldDouble(ON_NODES, ONE_TIME)
+        srcField.setNature(IntensiveMaximum)
+        srcField.setMesh(s); srcField.setName("field")
+        srcField.setArray(DataArrayDouble([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]))
+        tgtF = aRemapper.transferField(srcField, 1e+300)
+        srcField.writeVTK("/tmp/src3d.vtu")
+        tgtF.writeVTK("/tmp/tgt3d.vtu")
+        ref = []
+        val = tgtF.getArray().getValues()
+        for i, ref_v in enumerate(ref):
+            self.assertAlmostEqual(ref_v, val[i])        
+        pass
+
     @unittest.skipUnless(MEDCouplingHasNumPyBindings() and MEDCouplingHasSciPyBindings(),"requires numpy AND scipy")
     def testGetCrudeCSRMatrix1(self):
         """ testing CSR matrix output using numpy/scipy.
@@ -945,4 +1052,5 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         pass
     pass
 
-unittest.main()
+if __name__ == "__main__":
+  unittest.main()