*** empty log message ***

diff --git a/src/INTERP_KERNEL/IntegralUniformIntersector.hxx b/src/INTERP_KERNEL/IntegralUniformIntersector.hxx
new file mode 100644 (file)
index 0000000..8297323
--- /dev/null
+++ b/src/INTERP_KERNEL/IntegralUniformIntersector.hxx
@@ -0,0 +1,70 @@
+//  Copyright (C) 2007-2008  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 __INTEGRALUNIFORMINTERSECTOR_HXX__
+#define __INTEGRALUNIFORMINTERSECTOR_HXX__
+
+#include "TargetIntersector.hxx"
+
+#include <cmath>
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  class IntegralUniformIntersector : public TargetIntersector<MyMeshType,MyMatrix>
+  {
+  public:
+    typedef typename MyMeshType::MyConnType ConnType;
+  public:
+    IntegralUniformIntersector(const MyMeshType& mesh, bool isAbs);
+    double performNormalization(double val) const { if(_is_abs) return fabs(val); else return val; }
+    void setFromTo(bool val) { _from_to=val; }
+    void putValueIn(ConnType i, double val, MyMatrix& res) const;
+  protected:
+    const MyMeshType& _mesh;
+    //! if false means fromIntegralUniform if true means toIntegralUniform
+    bool _from_to;
+    bool _is_abs;
+  };
+
+  template<class MyMeshType, class MyMatrix>
+  class IntegralUniformIntersectorP0 : public IntegralUniformIntersector<MyMeshType,MyMatrix>
+  {
+  public:
+    typedef typename MyMeshType::MyConnType ConnType;
+  public:
+    IntegralUniformIntersectorP0(const MyMeshType& mesh, bool isAbs);
+    int getNumberOfRowsOfResMatrix() const;
+    int getNumberOfColsOfResMatrix() const;
+    void intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res);
+  };
+
+  template<class MyMeshType, class MyMatrix>
+  class IntegralUniformIntersectorP1 : public IntegralUniformIntersector<MyMeshType,MyMatrix>
+  {
+  public:
+    typedef typename MyMeshType::MyConnType ConnType;
+  public:
+    IntegralUniformIntersectorP1(const MyMeshType& mesh, bool isAbs);
+    int getNumberOfRowsOfResMatrix() const;
+    int getNumberOfColsOfResMatrix() const;
+    void intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res);
+  };
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/IntegralUniformIntersector.txx b/src/INTERP_KERNEL/IntegralUniformIntersector.txx
new file mode 100644 (file)
index 0000000..0acc2a1
--- /dev/null
+++ b/src/INTERP_KERNEL/IntegralUniformIntersector.txx
@@ -0,0 +1,155 @@
+//  Copyright (C) 2007-2008  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 __INTEGRALUNIFORMINTERSECTOR_TXX__
+#define __INTEGRALUNIFORMINTERSECTOR_TXX__
+
+#include "IntegralUniformIntersector.hxx"
+#include "VolSurfUser.txx"
+
+namespace INTERP_KERNEL
+{
+  template<class MyMeshType, class MyMatrix>
+  IntegralUniformIntersector<MyMeshType,MyMatrix>::IntegralUniformIntersector(const MyMeshType& mesh, bool isAbs):_mesh(mesh),_from_to(false),_is_abs(isAbs)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  void IntegralUniformIntersector<MyMeshType,MyMatrix>::putValueIn(ConnType iInCMode, double val1, MyMatrix& res) const
+  {
+    static const NumberingPolicy numPol=MyMeshType::My_numPol;
+    double val=performNormalization(val1);
+    if(_from_to)
+      {
+        typename MyMatrix::value_type& resRow=res[0];
+        typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(iInCMode));
+        if(iterRes==resRow.end())
+          resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(iInCMode),val));
+        else
+          {
+            double val2=(*iterRes).second+val;
+            resRow.erase(OTT<ConnType,numPol>::indFC(iInCMode));
+            resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(iInCMode),val2));
+          }
+      }
+    else
+      {
+        typename MyMatrix::value_type& resRow=res[iInCMode];
+        typename MyMatrix::value_type::const_iterator iterRes=resRow.find(OTT<ConnType,numPol>::indFC(0));
+        if(iterRes==resRow.end())
+          resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(0),val));
+        else
+          {
+            double val2=(*iterRes).second+val;
+            resRow.erase(OTT<ConnType,numPol>::indFC(0));
+            resRow.insert(std::make_pair(OTT<ConnType,numPol>::indFC(0),val2));
+          }
+      }
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  IntegralUniformIntersectorP0<MyMeshType,MyMatrix>::IntegralUniformIntersectorP0(const MyMeshType& mesh, bool isAbs):IntegralUniformIntersector<MyMeshType,MyMatrix>(mesh,isAbs)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int IntegralUniformIntersectorP0<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
+  {
+    if(IntegralUniformIntersector<MyMeshType,MyMatrix>::_from_to)
+      return 1;
+    else
+      return IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfElements();
+  }
+  
+  template<class MyMeshType, class MyMatrix>
+  int IntegralUniformIntersectorP0<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
+  {
+    if(IntegralUniformIntersector<MyMeshType,MyMatrix>::_from_to)
+      return IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfElements();
+    else
+      return 1;
+  }
+  
+  template<class MyMeshType, class MyMatrix>
+  void IntegralUniformIntersectorP0<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
+  {
+    static const NumberingPolicy numPol=MyMeshType::My_numPol;
+    res.resize(getNumberOfRowsOfResMatrix());
+    unsigned long nbelem=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfElements();
+    const ConnType *connIndx=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getConnectivityIndexPtr();
+    const ConnType *conn=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getConnectivityPtr();
+    const double *coords=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getCoordinatesPtr();
+    for(unsigned long i=0;i<nbelem;i++)
+      {
+        INTERP_KERNEL::NormalizedCellType t=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
+        double val=computeVolSurfOfCell<ConnType,numPol,MyMeshType::MY_SPACEDIM>(t,conn+OTT<ConnType,numPol>::ind2C(connIndx[i]),connIndx[i+1]-connIndx[i],coords);
+        IntegralUniformIntersector<MyMeshType,MyMatrix>::putValueIn(i,val,res);
+      }
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  IntegralUniformIntersectorP1<MyMeshType,MyMatrix>::IntegralUniformIntersectorP1(const MyMeshType& mesh, bool isAbs):IntegralUniformIntersector<MyMeshType,MyMatrix>(mesh,isAbs)
+  {
+  }
+
+  template<class MyMeshType, class MyMatrix>
+  int IntegralUniformIntersectorP1<MyMeshType,MyMatrix>::getNumberOfRowsOfResMatrix() const
+  {
+    if(IntegralUniformIntersector<MyMeshType,MyMatrix>::_from_to)
+      return 1;
+    else
+      return IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfNodes();
+  }
+  
+  template<class MyMeshType, class MyMatrix>
+  int IntegralUniformIntersectorP1<MyMeshType,MyMatrix>::getNumberOfColsOfResMatrix() const
+  {
+    if(IntegralUniformIntersector<MyMeshType,MyMatrix>::_from_to)
+      return IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfNodes();
+    else
+      return 1;
+  }
+  
+  template<class MyMeshType, class MyMatrix>
+  void IntegralUniformIntersectorP1<MyMeshType,MyMatrix>::intersectCells(ConnType targetCell, const std::vector<ConnType>& srcCells, MyMatrix& res)
+  {
+    static const NumberingPolicy numPol=MyMeshType::My_numPol;
+    res.resize(getNumberOfRowsOfResMatrix());
+    unsigned long nbelem=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getNumberOfElements();
+    const ConnType *connIndx=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getConnectivityIndexPtr();
+    const ConnType *conn=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getConnectivityPtr();
+    const double *coords=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getCoordinatesPtr();
+    for(unsigned long i=0;i<nbelem;i++)
+      {
+        INTERP_KERNEL::NormalizedCellType t=IntegralUniformIntersector<MyMeshType,MyMatrix>::_mesh.getTypeOfElement(OTT<ConnType,numPol>::indFC(i));
+        int lgth=connIndx[i+1]-connIndx[i];
+        const ConnType *locConn=conn+OTT<ConnType,numPol>::ind2C(connIndx[i]);
+        double val=computeVolSurfOfCell<ConnType,numPol,MyMeshType::MY_SPACEDIM>(t,locConn,lgth,coords);
+        if(t==NORM_TRI3)
+          val/=3.;
+        else if(t==NORM_TETRA4)
+          val/=4.;
+        else
+          throw INTERP_KERNEL::Exception("Invalid cell type detected : must be TRI3 or TETRA4 ! ");
+        for(int j=0;j<lgth;j++)
+          IntegralUniformIntersector<MyMeshType,MyMatrix>::putValueIn(OTT<ConnType,numPol>::coo2C(locConn[j]),val,res);
+      }
+  }
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/Interpolation.hxx b/src/INTERP_KERNEL/Interpolation.hxx
index 57edda88e39dd51910f80ec76459af8a566057e0..fe34d1b028288b49709b39e815632d74152d5f96 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation.hxx
+++ b/src/INTERP_KERNEL/Interpolation.hxx
@@ -38,8 +38,15 @@ namespace INTERP_KERNEL
     Interpolation(const InterpolationOptions& io) :InterpolationOptions(io){}
     //interpolation of two triangular meshes.
     template<class MatrixType, class MyMeshType>
-    int interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result)
-    { return asLeaf().interpolateMeshes(mesh1,mesh2,result); }
+    int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result)
+    { return asLeaf().interpolateMeshes(meshS,meshT,result); }
+    template<class MyMeshType, class MatrixType>
+    int fromIntegralUniform(const MyMeshType& meshT, MatrixType& result, const char *method) { return fromToIntegralUniform(true,meshT,result,method); }
+    template<class MyMeshType, class MatrixType>
+    int toIntegralUniform(const MyMeshType& meshS, MatrixType& result, const char *method) { return fromToIntegralUniform(false,meshS,result,method); }
+  protected:
+    template<class MyMeshType, class MatrixType>
+    int fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method);
   protected:
     TrueMainInterpolator& asLeaf() { return static_cast<TrueMainInterpolator&>(*this); }
   };

diff --git a/src/INTERP_KERNEL/Interpolation.txx b/src/INTERP_KERNEL/Interpolation.txx
new file mode 100644 (file)
index 0000000..483d8d1
--- /dev/null
+++ b/src/INTERP_KERNEL/Interpolation.txx
@@ -0,0 +1,58 @@
+//  Copyright (C) 2007-2008  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 __INTERPOLATION_TXX__
+#define __INTERPOLATION_TXX__
+
+#include "Interpolation.hxx"
+#include "IntegralUniformIntersector.hxx"
+#include "IntegralUniformIntersector.txx"
+
+namespace INTERP_KERNEL
+{ 
+  template<class TrueMainInterpolator>
+  template<class MyMeshType, class MatrixType>
+  int Interpolation<TrueMainInterpolator>::fromToIntegralUniform(bool fromTo, const MyMeshType& mesh, MatrixType& result, const char *method)
+  {
+    typedef typename MyMeshType::MyConnType ConnType;
+    std::string methodCPP(method);
+    int ret=-1;
+    if(methodCPP=="P0")
+      {
+        IntegralUniformIntersectorP0<MyMeshType,MatrixType> intersector(mesh,InterpolationOptions::getMeasureAbsStatus());
+        intersector.setFromTo(fromTo);
+        std::vector<ConnType> tmp;
+        intersector.intersectCells(0,tmp,result);
+        ret=intersector.getNumberOfColsOfResMatrix();
+      }
+    else if(methodCPP=="P1")
+      {
+        IntegralUniformIntersectorP1<MyMeshType,MatrixType> intersector(mesh,InterpolationOptions::getMeasureAbsStatus());
+        intersector.setFromTo(fromTo);
+        std::vector<ConnType> tmp;
+        intersector.intersectCells(0,tmp,result);
+        ret=intersector.getNumberOfColsOfResMatrix();
+      }
+    else
+      throw INTERP_KERNEL::Exception("Invalid method specified in fromIntegralUniform : must be in { \"P0\", \"P1\"}");
+    return ret;
+  }
+}
+
+#endif
+

diff --git a/src/INTERP_KERNEL/Interpolation3D.txx b/src/INTERP_KERNEL/Interpolation3D.txx
index d24cbdc8adf36b582d2d4de8fe4b4a3003ee2e19..35ee63b3417a9921082eed72192a817be03b29b7 100644 (file)
--- a/src/INTERP_KERNEL/Interpolation3D.txx
+++ b/src/INTERP_KERNEL/Interpolation3D.txx
@@ -20,6 +20,7 @@
 #define __INTERPOLATION3D_TXX__
 
 #include "Interpolation3D.hxx"
+#include "Interpolation.txx"
 #include "MeshElement.txx"
 #include "TransformedTriangle.hxx"
 #include "PolyhedronIntersector.txx"
@@ -315,7 +316,6 @@ namespace INTERP_KERNEL
     return ret;
 
   }
-
 }
 
 #endif

diff --git a/src/INTERP_KERNEL/InterpolationOptions.cxx b/src/INTERP_KERNEL/InterpolationOptions.cxx
index 6bbb520b60850e318f36d6db21777876391a400d..c4514b16009ec816a4b3c32b1286d5c4040703af 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationOptions.cxx
+++ b/src/INTERP_KERNEL/InterpolationOptions.cxx
@@ -40,6 +40,8 @@ const char INTERP_KERNEL::InterpolationOptions::DO_ROTATE_STR[]="DoRotate";
 
 const char INTERP_KERNEL::InterpolationOptions::ORIENTATION_STR[]="Orientation";
 
+const char INTERP_KERNEL::InterpolationOptions::MEASURE_ABS_STR[]="MeasureAbs";
+
 const char INTERP_KERNEL::InterpolationOptions::INTERSEC_TYPE_STR[]="IntersectionType";
 
 const char INTERP_KERNEL::InterpolationOptions::SPLITTING_POLICY_STR[]="SplittingPolicy";
@@ -106,6 +108,11 @@ bool INTERP_KERNEL::InterpolationOptions::setOptionInt(const std::string& key, i
         setOrientation(value);
         return true;
       }
+    else if(key==MEASURE_ABS_STR)
+      {
+        bool valBool=(value!=0);
+        setMeasureAbsStatus(valBool);
+      }
     else
       return false;
 }

diff --git a/src/INTERP_KERNEL/InterpolationOptions.hxx b/src/INTERP_KERNEL/InterpolationOptions.hxx
index bfc8435080857c365e206a9833b616807b6794f5..89430d1f0d871613c9898211e94e1eed0c928681 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationOptions.hxx
+++ b/src/INTERP_KERNEL/InterpolationOptions.hxx
@@ -46,6 +46,7 @@ namespace INTERP_KERNEL
     double _bounding_box_adjustment_abs ;
     double _max_distance_for_3Dsurf_intersect;
     int _orientation ;
+    bool _measure_abs;
     SplittingPolicy _splitting_policy ;
 
   public:
@@ -76,6 +77,9 @@ namespace INTERP_KERNEL
 
     int getOrientation() const { return _orientation; }
     void setOrientation(int o) { _orientation=o; }
+
+    bool getMeasureAbsStatus() const { return _measure_abs; }
+    void setMeasureAbsStatus(bool newStatus) { _measure_abs=newStatus; }
     
     SplittingPolicy getSplittingPolicy() const { return _splitting_policy; }
     void setSplittingPolicy(SplittingPolicy sp) { _splitting_policy=sp; }
@@ -90,6 +94,7 @@ namespace INTERP_KERNEL
       _bounding_box_adjustment_abs=0.;
       _max_distance_for_3Dsurf_intersect=DFT_MAX_DIST_3DSURF_INTERSECT;
       _orientation=0;
+      _measure_abs=true;
       _splitting_policy=GENERAL_48;
     }
     bool setOptionDouble(const std::string& key, double value);
@@ -108,6 +113,7 @@ namespace INTERP_KERNEL
     static const char PRINT_LEV_STR[];
     static const char DO_ROTATE_STR[];
     static const char ORIENTATION_STR[];
+    static const char MEASURE_ABS_STR[];
     static const char INTERSEC_TYPE_STR[];
     static const char SPLITTING_POLICY_STR[];
     static const char TRIANGULATION_INTERSECT2D_STR[];

diff --git a/src/INTERP_KERNEL/InterpolationPlanar.hxx b/src/INTERP_KERNEL/InterpolationPlanar.hxx
index f4321c95ed74af39bd2abf8404fec0c70ca2a911..51bd627734dce423628d709ccf13846ceb6dd6d7 100755 (executable)
--- a/src/INTERP_KERNEL/InterpolationPlanar.hxx
+++ b/src/INTERP_KERNEL/InterpolationPlanar.hxx
@@ -31,8 +31,6 @@ namespace INTERP_KERNEL
   {
   private:
     double _dim_caracteristic;
-    static const double DEFAULT_PRECISION;
-
   public:
     InterpolationPlanar();
     InterpolationPlanar(const InterpolationOptions & io);
@@ -43,15 +41,13 @@ namespace INTERP_KERNEL
     
     // Main function to interpolate triangular and quadratic meshes
     template<class MyMeshType, class MatrixType>
-      int interpolateMeshes(const MyMeshType& mesh1, const MyMeshType& mesh2, MatrixType& result, const char *method);
-    
+    int interpolateMeshes(const MyMeshType& meshS, const MyMeshType& meshT, MatrixType& result, const char *method);
   public:
     bool doRotate() const { return asLeafInterpPlanar().doRotate(); }
     double medianPlane() const { return asLeafInterpPlanar().medianPlane(); }
     template<class MyMeshType, class MyMatrixRow>
       void performAdjustmentOfBB(PlanarIntersector<MyMeshType,MyMatrixRow>* intersector, std::vector<double>& bbox) const
     { return asLeafInterpPlanar().performAdjustmentOfBB(intersector,bbox); }
-
   protected:
     RealPlanar& asLeafInterpPlanar() { return static_cast<RealPlanar&>(*this); }
     const RealPlanar& asLeafInterpPlanar() const { return static_cast< const RealPlanar& >(*this); }

diff --git a/src/INTERP_KERNEL/InterpolationPlanar.txx b/src/INTERP_KERNEL/InterpolationPlanar.txx
index 1e720fd25fdf2b451e9badb0c37ecdd7ae24264b..49eed13aad29a391144e1fda37ab49b953a59aed 100644 (file)
--- a/src/INTERP_KERNEL/InterpolationPlanar.txx
+++ b/src/INTERP_KERNEL/InterpolationPlanar.txx
@@ -20,6 +20,7 @@
 #define __INTERPOLATIONPLANAR_TXX__
 
 #include "InterpolationPlanar.hxx"
+#include "Interpolation.txx"
 #include "InterpolationOptions.hxx"
 #include "PlanarIntersector.hxx"
 #include "PlanarIntersector.txx"
@@ -36,10 +37,6 @@
 
 namespace INTERP_KERNEL
 {
-
-  template<class RealPlanar>
-  const double InterpolationPlanar<RealPlanar>::DEFAULT_PRECISION=1.e-12;
-
   /**
    * \defgroup interpolationPlanar InterpolationPlanar
    *
@@ -59,7 +56,6 @@ namespace INTERP_KERNEL
   {
   }
 
-
   /**
    *  \brief  Function used to set the options for the intersection calculation
    * \details The following options can be modified:

diff --git a/src/INTERP_KERNEL/Makefile.am b/src/INTERP_KERNEL/Makefile.am
index 0403e79f9422e96976bea0095b1c3be9a7d97e1f..042c2c49afff3afc4eb36073cbc06bcc2de244e1 100644 (file)
--- a/src/INTERP_KERNEL/Makefile.am
+++ b/src/INTERP_KERNEL/Makefile.am
@@ -41,6 +41,7 @@ Geometric2DIntersector.txx            \
 INTERPKERNELDefines.hxx                        \
 InterpKernelMatrix.hxx                 \
 Interpolation.hxx                      \
+Interpolation.txx                       \
 Interpolation2D.hxx                    \
 Interpolation2D.txx                    \
 Interpolation3D.hxx                    \
@@ -91,6 +92,8 @@ TransformedTriangleInline.hxx         \
 TranslationRotationMatrix.hxx          \
 TriangulationIntersector.hxx           \
 TriangulationIntersector.txx           \
+IntegralUniformIntersector.hxx          \
+IntegralUniformIntersector.txx          \
 UnitTetraIntersectionBary.hxx          \
 VTKNormalizedUnstructuredMesh.hxx      \
 VTKNormalizedUnstructuredMesh.txx      \
@@ -104,7 +107,9 @@ PlanarIntersectorP0P1.hxx   \
 PlanarIntersectorP0P1.txx      \
 PlanarIntersectorP1P0.hxx      \
 PlanarIntersectorP1P0.txx      \
-VolSurfFormulae.hxx
+VolSurfFormulae.hxx             \
+VolSurfUser.hxx                 \
+VolSurfUser.txx
 
 # Libraries targets
 

diff --git a/src/INTERP_KERNEL/VolSurfFormulae.hxx b/src/INTERP_KERNEL/VolSurfFormulae.hxx
index 52c49856ab635b8a795fcd896fe5e6f76228a023..7effb8eff760bd3dc16abcaa6731a0894b0aad5b 100644 (file)
--- a/src/INTERP_KERNEL/VolSurfFormulae.hxx
+++ b/src/INTERP_KERNEL/VolSurfFormulae.hxx
@@ -16,8 +16,8 @@
 //
 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-#ifndef VOLSURFFORMULAE
-#define VOLSURFFORMULAE
+#ifndef __VOLSURFFORMULAE_HXX__
+#define __VOLSURFFORMULAE_HXX__
 
 #include <math.h>
 

diff --git a/src/INTERP_KERNEL/VolSurfUser.hxx b/src/INTERP_KERNEL/VolSurfUser.hxx
new file mode 100644 (file)
index 0000000..ef89fcd
--- /dev/null
+++ b/src/INTERP_KERNEL/VolSurfUser.hxx
@@ -0,0 +1,33 @@
+//  Copyright (C) 2007-2008  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 __VOLSURFUSER_HXX__
+#define __VOLSURFUSER_HXX__
+
+#include "NormalizedUnstructuredMesh.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class ConnType, NumberingPolicy numPolConn, int SPACEDIM>
+  double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords);
+
+  template<class ConnType, NumberingPolicy numPolConn>
+  double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim);
+}
+
+#endif

diff --git a/src/INTERP_KERNEL/VolSurfUser.txx b/src/INTERP_KERNEL/VolSurfUser.txx
new file mode 100644 (file)
index 0000000..e657796
--- /dev/null
+++ b/src/INTERP_KERNEL/VolSurfUser.txx
@@ -0,0 +1,168 @@
+//  Copyright (C) 2007-2008  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 __VOLSURFUSER_TXX__
+#define __VOLSURFUSER_TXX__
+
+#include "VolSurfUser.hxx"
+#include "VolSurfFormulae.hxx"
+#include "InterpolationUtils.hxx"
+
+namespace INTERP_KERNEL
+{
+  template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
+  double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords)
+  {
+    switch(type)
+      {
+      case INTERP_KERNEL::NORM_TRI3 :
+      case INTERP_KERNEL::NORM_TRI6 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+              
+          return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
+                                                     coords+(SPACEDIM*N2),
+                                                     coords+(SPACEDIM*N3),
+                                                     SPACEDIM);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_QUAD4 :
+      case INTERP_KERNEL::NORM_QUAD8 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+              
+          return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
+                                                     coords+SPACEDIM*N2,
+                                                     coords+SPACEDIM*N3,
+                                                     coords+SPACEDIM*N4,
+                                                     SPACEDIM);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_POLYGON :
+        {          
+          const double **pts=new const double *[lgth];
+          for(int inod=0;inod<lgth;inod++)
+            pts[inod] = coords+3*OTT<ConnType,numPol>::coo2C(connec[inod]);
+          double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM);
+          delete [] pts;
+          return val;
+        }
+        break;
+      case INTERP_KERNEL::NORM_TETRA4 :
+      case INTERP_KERNEL::NORM_TETRA10 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+              
+          return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
+                                                        coords+SPACEDIM*N2,
+                                                        coords+SPACEDIM*N3,
+                                                        coords+SPACEDIM*N4);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_PYRA5 :
+      case INTERP_KERNEL::NORM_PYRA13 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+              
+          return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
+                                                       coords+SPACEDIM*N2,
+                                                       coords+SPACEDIM*N3,
+                                                       coords+SPACEDIM*N4,
+                                                       coords+SPACEDIM*N5);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_PENTA6 :
+      case INTERP_KERNEL::NORM_PENTA15 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+          int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
+              
+          return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
+                                                        coords+SPACEDIM*N2,
+                                                        coords+SPACEDIM*N3,
+                                                        coords+SPACEDIM*N4,
+                                                        coords+SPACEDIM*N5,
+                                                        coords+SPACEDIM*N6);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_HEXA8 :
+      case INTERP_KERNEL::NORM_HEXA20 :
+        {
+          int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
+          int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
+          int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
+          int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
+          int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
+          int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
+          int N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
+          int N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
+              
+          return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
+                                                       coords+SPACEDIM*N2,
+                                                       coords+SPACEDIM*N3,
+                                                       coords+SPACEDIM*N4,
+                                                       coords+SPACEDIM*N5,
+                                                       coords+SPACEDIM*N6,
+                                                       coords+SPACEDIM*N7,
+                                                       coords+SPACEDIM*N8);
+        }
+        break;
+            
+      case INTERP_KERNEL::NORM_POLYHED :
+        {
+          throw INTERP_KERNEL::Exception("Polyedra Not yet implemented !");
+        }
+        break;
+      default:
+        throw INTERP_KERNEL::Exception("Not recognized cell type to get Area/Volume on it !");
+      }
+  }
+
+  template<class ConnType, NumberingPolicy numPolConn>
+  double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
+  {
+    if(spaceDim==3)
+      return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
+    if(spaceDim==2)
+      return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
+    throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 2 or 3");
+  }
+}
+
+#endif

diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx
index 33f7214c05bea017636a70720ec110b40ac7241d..19deae8d151ea169cb295d951594f827332cb829 100644 (file)
--- a/src/MEDCoupling/MEDCouplingCMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.cxx
@@ -175,7 +175,7 @@ void MEDCouplingCMesh::getBoundingBox(double *bbox) const
   //not implemented yet !
 }
 
-MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField() const
+MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 {
   //not implemented yet !
   return 0;

diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx
index 212535c441219b2981dfe928156cfe8274e4cdf1..59fc9c2e7ca62eb4405444edc7a296e61396e16d 100644 (file)
--- a/src/MEDCoupling/MEDCouplingCMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.hxx
@@ -45,7 +45,7 @@ namespace ParaMEDMEM
                    DataArrayDouble *coordsZ=0);
     // tools
     void getBoundingBox(double *bbox) const;
-    MEDCouplingFieldDouble *getMeasureField() const;
+    MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
   private:
     MEDCouplingCMesh();
     ~MEDCouplingCMesh();

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 5b3e5cac6106ea7f7c568ccfa1a8ab50fb6c8de3..d1c82fc26dff59d13d98d2716538e930989f7fea 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -86,9 +86,9 @@ void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMe
     }
 }
 
-MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh) const
+MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const
 {
-  return mesh->getMeasureField();
+  return mesh->getMeasureField(isAbs);
 }
 
 /*!
@@ -149,7 +149,7 @@ void MEDCouplingFieldDiscretizationP1::checkCoherencyBetween(const MEDCouplingMe
     }
 }
 
-MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh) const
+MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const
 {
   //not implemented yet.
   //Dual mesh to build

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index 6061bf20f883e2af6a0193416935cf9d10eb73ff..f035ba3edfe284d9570d0d1e4d7bed418075ad4c 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -42,7 +42,7 @@ namespace ParaMEDMEM
     virtual int getNumberOfTuples(const MEDCouplingMesh *mesh) const = 0;
     virtual void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception) = 0;
     virtual void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception) = 0;
-    virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const = 0;
+    virtual MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const = 0;
     virtual MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const = 0;
   };
 
@@ -56,7 +56,7 @@ namespace ParaMEDMEM
     int getNumberOfTuples(const MEDCouplingMesh *mesh) const;
     void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
     void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
-    MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+    MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
     MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
   public:
     static const char REPR[];
@@ -73,7 +73,7 @@ namespace ParaMEDMEM
     int getNumberOfTuples(const MEDCouplingMesh *mesh) const;
     void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
     void checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception);
-    MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh) const;
+    MEDCouplingFieldDouble *getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const;
     MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, const MEDCouplingMesh *mesh, DataArrayInt *&di) const;
     static DataArrayInt *invertArrayO2N2N2O(const MEDCouplingMesh *mesh, const DataArrayInt *di);
   public:

diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
index 6b1fa4e64e3dbf46034620806600d74bd3b388f8..5c673c463fadbfef7ceca9ca06280b33552c5b75 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
@@ -89,11 +89,11 @@ double MEDCouplingFieldDouble::accumulate(int compId) const
   return ret;
 }
 
-double MEDCouplingFieldDouble::measureAccumulate(int compId) const
+double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const
 {
   if(!_mesh)
     throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate");
-  MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh);
+  MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh,isWAbs);
   const double *ptr=weight->getArray()->getConstPointer();
   int nbOfValues=weight->getArray()->getNbOfElems();
   double ret=0.;

diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx
index a05035b0f3be75ca46456dfc8d96c98a264a9488..0f37543affbcabcfc8cae1dc57467a0bd0d3676f 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx
@@ -43,7 +43,7 @@ namespace ParaMEDMEM
     void setArray(DataArrayDouble *array);
     DataArrayDouble *getArray() const { return _time_discr->getArray(); }
     double accumulate(int compId) const;
-    double measureAccumulate(int compId) const;
+    double measureAccumulate(int compId, bool isWAbs) const;
     void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception);
     void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception);
     //! \b temporary

diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx
index 61324c2f8b15c7d02b6ad9e89d8b7c3caa7b3952..597088efe39fb96a27006bd42aef2d66972e04a5 100644 (file)
--- a/src/MEDCoupling/MEDCouplingMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingMesh.hxx
@@ -50,7 +50,7 @@ namespace ParaMEDMEM
     virtual int getMeshDimension() const = 0;
     // tools
     virtual void getBoundingBox(double *bbox) const = 0;
-    virtual MEDCouplingFieldDouble *getMeasureField() const = 0;
+    virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0;
   protected:
     MEDCouplingMesh() { }
     MEDCouplingMesh(const MEDCouplingMesh& other):_name(other._name) { }

diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 4fa962cb2dba8fa24b9f986e0bafa22b67153a9c..a560bef915d106a072b8d39562c31341c45ae28d 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -19,7 +19,7 @@
 #include "MEDCouplingUMesh.hxx"
 #include "MEDCouplingFieldDouble.hxx"
 #include "CellModel.hxx"
-#include "VolSurfFormulae.hxx"
+#include "VolSurfUser.txx"
 
 #include <sstream>
 #include <limits>
@@ -755,185 +755,29 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfKeepCoords(const int *start
  * param field field on which cells the volumes are required
  * return field containing the volumes, area or length depending the meshdimension.
  */
-MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField() const
-{
-  int ipt, type ;
-  int nbelem       = getNumberOfCells() ;
-  int dim_mesh     = getMeshDimension();
-  int dim_space    = getSpaceDimension() ;
-  const double *coords = getCoords()->getConstPointer() ;
-  const int *connec = getNodalConnectivity()->getConstPointer() ;
-  const int *connec_index = getNodalConnectivityIndex()->getConstPointer() ;
-
-
-  MEDCouplingFieldDouble* field = MEDCouplingFieldDouble::New(ON_CELLS);
-  DataArrayDouble* array = DataArrayDouble::New() ;
-  array->alloc(nbelem, 1) ;
-  double *area_vol = array->getPointer() ;
-
-  switch (dim_mesh)
+MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const
+{
+  int ipt;
+  INTERP_KERNEL::NormalizedCellType type;
+  int nbelem=getNumberOfCells();
+  int dim_space=getSpaceDimension();
+  const double *coords=getCoords()->getConstPointer();
+  const int *connec=getNodalConnectivity()->getConstPointer();
+  const int *connec_index=getNodalConnectivityIndex()->getConstPointer();
+
+  MEDCouplingFieldDouble* field=MEDCouplingFieldDouble::New(ON_CELLS);
+  DataArrayDouble* array=DataArrayDouble::New();
+  array->alloc(nbelem,1);
+  double *area_vol = array->getPointer();
+  for(int iel=0;iel<nbelem;iel++)
     {
-    case 2: // getting the areas
-      for ( int iel=0 ; iel<nbelem ; iel++ )
-        {
-          ipt = connec_index[iel] ;
-          type = connec[ipt] ;
-
-          switch ( type )
-            {
-            case INTERP_KERNEL::NORM_TRI3 :
-            case INTERP_KERNEL::NORM_TRI6 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-
-                area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForTria(coords+(dim_space*N1),
-                                                                       coords+(dim_space*N2),
-                                                                       coords+(dim_space*N3),
-                                                                       dim_space));
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_QUAD4 :
-            case INTERP_KERNEL::NORM_QUAD8 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-                int N4 = connec[ipt+4];
-
-                area_vol[iel]=fabs(INTERP_KERNEL::calculateAreaForQuad(coords+dim_space*N1,
-                                                                       coords+dim_space*N2,
-                                                                       coords+dim_space*N3,
-                                                                       coords+dim_space*N4,
-                                                                       dim_space));
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_POLYGON :
-              {
-                // We must remember that the first item is the type. That's
-                // why we substract 1 to get the number of nodes of this polygon
-                int size = connec_index[iel+1] - connec_index[iel] - 1 ;
-
-                const double **pts = new const double *[size] ;
-
-                for ( int inod=0 ; inod<size ; inod++ )
-                  {
-                    // Remember the first item is the type
-                    pts[inod] = coords+dim_space*connec[ipt+inod+1] ;
-                  }
-
-                area_vol[iel]=INTERP_KERNEL::calculateAreaForPolyg(pts,size,dim_space);
-                delete [] pts;
-              }
-              break ;
-
-            default :
-              throw INTERP_KERNEL::Exception("Bad Support to get Areas on it !");
-
-            } // End of switch
-
-        } // End of the loop over the cells
-      break;
-    case 3: // getting the volumes
-      for ( int iel=0 ; iel<nbelem ; iel++ )
-        {
-          ipt = connec_index[iel] ;
-          type = connec[ipt] ;
-
-          switch ( type )
-            {
-            case INTERP_KERNEL::NORM_TETRA4 :
-            case INTERP_KERNEL::NORM_TETRA10 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-                int N4 = connec[ipt+4];
-
-                area_vol[iel]=INTERP_KERNEL::calculateVolumeForTetra(coords+dim_space*N1,
-                                                                     coords+dim_space*N2,
-                                                                     coords+dim_space*N3,
-                                                                     coords+dim_space*N4) ;
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_PYRA5 :
-            case INTERP_KERNEL::NORM_PYRA13 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-                int N4 = connec[ipt+4];
-                int N5 = connec[ipt+5];
-
-                area_vol[iel]=INTERP_KERNEL::calculateVolumeForPyra(coords+dim_space*N1,
-                                                                    coords+dim_space*N2,
-                                                                    coords+dim_space*N3,
-                                                                    coords+dim_space*N4,
-                                                                    coords+dim_space*N5) ;
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_PENTA6 :
-            case INTERP_KERNEL::NORM_PENTA15 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-                int N4 = connec[ipt+4];
-                int N5 = connec[ipt+5];
-                int N6 = connec[ipt+6];
-
-                area_vol[iel]=INTERP_KERNEL::calculateVolumeForPenta(coords+dim_space*N1,
-                                                                     coords+dim_space*N2,
-                                                                     coords+dim_space*N3,
-                                                                     coords+dim_space*N4,
-                                                                     coords+dim_space*N5,
-                                                                     coords+dim_space*N6) ;
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_HEXA8 :
-            case INTERP_KERNEL::NORM_HEXA20 :
-              {
-                int N1 = connec[ipt+1];
-                int N2 = connec[ipt+2];
-                int N3 = connec[ipt+3];
-                int N4 = connec[ipt+4];
-                int N5 = connec[ipt+5];
-                int N6 = connec[ipt+6];
-                int N7 = connec[ipt+7];
-                int N8 = connec[ipt+8];
-
-                area_vol[iel]=INTERP_KERNEL::calculateVolumeForHexa(coords+dim_space*N1,
-                                                                    coords+dim_space*N2,
-                                                                    coords+dim_space*N3,
-                                                                    coords+dim_space*N4,
-                                                                    coords+dim_space*N5,
-                                                                    coords+dim_space*N6,
-                                                                    coords+dim_space*N7,
-                                                                    coords+dim_space*N8) ;
-              }
-              break ;
-
-            case INTERP_KERNEL::NORM_POLYHED :
-              {
-                throw INTERP_KERNEL::Exception("Not yet implemented !");
-              }
-              break ;
-
-            default:
-              throw INTERP_KERNEL::Exception("Bad Support to get Volume on it !");
-            }
-        }
-      break;
-    default:
-      throw INTERP_KERNEL::Exception("interpolation is not available for this dimension");
+      ipt=connec_index[iel];
+      type=(INTERP_KERNEL::NormalizedCellType)connec[ipt];
+      area_vol[iel]=INTERP_KERNEL::computeVolSurfOfCell2<int,INTERP_KERNEL::ALL_C_MODE>(type,connec+ipt+1,connec_index[iel+1]-ipt-1,coords,dim_space);
     }
-  
+  if(isAbs)
+    for(int iel=0;iel<nbelem;iel++)
+      area_vol[iel]=fabs(area_vol[iel]);
   field->setArray(array) ;
   array->decrRef();
   field->setMesh(const_cast<MEDCouplingUMesh *>(this));  

diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index 160a57ddbb0e7d726733fe458d7f368465b0e427..66000673d762b3b1f83759143d40d213bbebeb49 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -67,7 +67,7 @@ namespace ParaMEDMEM
     MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const;
     void renumberConnectivity(const int *newNodeNumbers);
     void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
-    MEDCouplingFieldDouble *getMeasureField() const;
+    MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
   private:
     MEDCouplingUMesh();
     MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy);

diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
index 19fa0abd5b0cf911df985b666ef5146786e9f269..fc20b2747aafe1d8eaba3f36dc2788ba9b1a50dd 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
@@ -242,7 +242,7 @@ void MEDCouplingUMeshDesc::renumberConnectivity(const int *newNodeNumbers)
 {
 }
 
-MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField() const
+MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField(bool isAbs) const
 {
   //not implemented yet.
   return 0;

diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
index b7e672435dc93400d382ceaaefc832ad86b5a61a..ca3c9b6d4087ebcba7c2ba33ec8fbdb484cd6433 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
@@ -53,7 +53,7 @@ namespace ParaMEDMEM
     void findBoundaryNodes(std::vector<int>& nodes) const;
     MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const;
     void renumberConnectivity(const int *newNodeNumbers);
-    MEDCouplingFieldDouble *getMeasureField() const;
+    MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
     DataArrayInt *zipCoordsTraducer();
   private:
     MEDCouplingUMeshDesc();

diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
index c0048929843e8af363cc66d4d185b85e659f3da2..ec673e640170868152ebb3da7149c458f57fc1ea 100644 (file)
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
@@ -1479,6 +1479,177 @@ void MEDCouplingBasicsTest::test3DInterpP0P0Empty()
   targetMesh->decrRef();
 }
 
+void MEDCouplingBasicsTest::test2DInterpP0IntegralUniform()
+{
+  MEDCouplingUMesh *targetMesh=build2DTargetMesh_1();
+  //
+  MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh);
+  INTERP_KERNEL::Interpolation2D myInterpolator;
+  vector<map<int,double> > res;
+  CPPUNIT_ASSERT_EQUAL(5,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+  res.clear();
+  CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[1][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[3][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[4][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+  res.clear();
+  targetMesh->decrRef();
+  //
+  targetMesh=build2DTargetMeshPerm_1();
+  MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper2(targetMesh);
+  INTERP_KERNEL::Interpolation2D myInterpolator2;
+  CPPUNIT_ASSERT(myInterpolator2.getMeasureAbsStatus());
+  CPPUNIT_ASSERT_EQUAL(5,myInterpolator2.fromIntegralUniform(targetWrapper2,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,sumAll(res),1e-12);
+  res.clear();
+  myInterpolator2.setMeasureAbsStatus(false);
+  CPPUNIT_ASSERT(!myInterpolator2.getMeasureAbsStatus());
+  CPPUNIT_ASSERT_EQUAL(5,myInterpolator2.fromIntegralUniform(targetWrapper2,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.125,res[0][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125,res[0][2],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25,res[0][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.75,sumAll(res),1e-12);
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DSurfInterpP0IntegralUniform()
+{
+  MEDCouplingUMesh *targetMesh=build3DSurfTargetMesh_1();
+  INTERP_KERNEL::Interpolation3DSurf myInterpolator;
+  MEDCouplingNormalizedUnstructuredMesh<3,2> targetWrapper(targetMesh);
+  vector<map<int,double> > res;
+  CPPUNIT_ASSERT_EQUAL(5,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[0][2],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][3],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][4],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12);
+  res.clear();
+  CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(5,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[1][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.125*sqrt(2.),res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[3][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.25*sqrt(2.),res[4][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.*sqrt(2.),sumAll(res),1e-12);
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DInterpP0IntegralUniform()
+{
+  MEDCouplingUMesh *targetMesh=build3DTargetMesh_1();
+  INTERP_KERNEL::Interpolation3D myInterpolator;
+  MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(targetMesh);
+  vector<map<int,double> > res;
+  CPPUNIT_ASSERT_EQUAL(8,myInterpolator.fromIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(125000.,res[0][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][1],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][2],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][3],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[0][4],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][5],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[0][6],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(3375000.,res[0][7],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6);
+  res.clear();
+  CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P0"));
+  CPPUNIT_ASSERT_EQUAL(8,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(125000.,res[0][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[1][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[2][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[3][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(375000.,res[4][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[5][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1125000.,res[6][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(3375000.,res[7][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6);
+  res.clear();
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test2DInterpP1IntegralUniform()
+{
+  MEDCouplingUMesh *targetMesh=build2DSourceMesh_1();
+  //
+  MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(targetMesh);
+  INTERP_KERNEL::Interpolation2D myInterpolator;
+  vector<map<int,double> > res;
+  CPPUNIT_ASSERT_EQUAL(4,myInterpolator.fromIntegralUniform(targetWrapper,res,"P1"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[0][1],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[0][2],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][3],1e-12);
+  res.clear();
+  CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P1"));
+  CPPUNIT_ASSERT_EQUAL(4,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[0][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[1][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.16666666666666666,res[2][0],1e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.33333333333333331,res[3][0],1e-12);
+  res.clear();
+  targetMesh->decrRef();
+}
+
+void MEDCouplingBasicsTest::test3DInterpP1IntegralUniform()
+{
+  MEDCouplingUMesh *sourceMesh=build3DSourceMesh_1();
+  //
+  MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(sourceMesh);
+  INTERP_KERNEL::Interpolation3D myInterpolator;
+  vector<map<int,double> > res;
+  CPPUNIT_ASSERT_EQUAL(9,myInterpolator.fromIntegralUniform(targetWrapper,res,"P1"));
+  CPPUNIT_ASSERT_EQUAL(1,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][1],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[0][2],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][3],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][4],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[0][5],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][6],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][7],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(2000000.,res[0][8],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6);
+  res.clear();
+  CPPUNIT_ASSERT_EQUAL(1,myInterpolator.toIntegralUniform(targetWrapper,res,"P1"));
+  CPPUNIT_ASSERT_EQUAL(9,(int)res.size());
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[0][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[1][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[2][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[3][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[4][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(500000.,res[5][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[6][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(833333.333333333,res[7][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(2000000.,res[8][0],1e-6);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(8000000.,sumAll(res),1e-6);
+  sourceMesh->decrRef();
+}
+
 MEDCouplingUMesh *MEDCouplingBasicsTest::build2DSourceMesh_1()
 {
   double sourceCoords[8]={-0.3,-0.3, 0.7,-0.3, -0.3,0.7, 0.7,0.7};
@@ -1517,6 +1688,27 @@ MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_1()
   return targetMesh;
 }
 
+MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMeshPerm_1()
+{
+  double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
+  int targetConn[18]={0,3,4,1, 1,2,4, 4,5,2, 6,7,4,3, 7,8,5,4};
+  MEDCouplingUMesh *targetMesh=MEDCouplingUMesh::New();
+  targetMesh->setMeshDimension(2);
+  targetMesh->allocateCells(5);
+  targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn);
+  targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+4);
+  targetMesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,targetConn+7);
+  targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+10);
+  targetMesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,targetConn+14);
+  targetMesh->finishInsertingCells();
+  DataArrayDouble *myCoords=DataArrayDouble::New();
+  myCoords->alloc(9,2);
+  std::copy(targetCoords,targetCoords+18,myCoords->getPointer());
+  targetMesh->setCoords(myCoords);
+  myCoords->decrRef();
+  return targetMesh;
+}
+
 MEDCouplingUMesh *MEDCouplingBasicsTest::build2DTargetMesh_2()
 {
   double targetCoords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };

diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index 49ddbc162e40c4af73bba9d1c1fb6c198facba5b..03941d92244b986b25c088e0475e3b0783fc739c 100644 (file)
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -63,6 +63,11 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( test3DInterpP0P1_1 );
     CPPUNIT_TEST( test3DInterpP1P0_1 );
     CPPUNIT_TEST( test3DInterpP0P0Empty );
+    CPPUNIT_TEST( test2DInterpP0IntegralUniform );
+    CPPUNIT_TEST( test3DSurfInterpP0IntegralUniform );
+    CPPUNIT_TEST( test3DInterpP0IntegralUniform );
+    CPPUNIT_TEST( test2DInterpP1IntegralUniform );
+    CPPUNIT_TEST( test3DInterpP1IntegralUniform );
     CPPUNIT_TEST_SUITE_END();
   public:
     void testArray();
@@ -97,9 +102,15 @@ namespace ParaMEDMEM
     void test3DInterpP0P1_1();
     void test3DInterpP1P0_1();
     void test3DInterpP0P0Empty();
+    void test2DInterpP0IntegralUniform();
+    void test3DSurfInterpP0IntegralUniform();
+    void test3DInterpP0IntegralUniform();
+    void test2DInterpP1IntegralUniform();
+    void test3DInterpP1IntegralUniform();
   private:
     MEDCouplingUMesh *build2DSourceMesh_1();
     MEDCouplingUMesh *build2DTargetMesh_1();
+    MEDCouplingUMesh *build2DTargetMeshPerm_1();
     MEDCouplingUMesh *build2DTargetMesh_2();
     MEDCouplingUMesh *build3DSurfSourceMesh_1();
     MEDCouplingUMesh *build3DSurfSourceMesh_2();

diff --git a/src/ParaMEDMEM/DEC.cxx b/src/ParaMEDMEM/DEC.cxx
index c564438bb6db7c2c4afdb385c0e7a8a011b878f2..8bb1f2d27fe3639d3d6f8058e9b82564c7a0397b 100644 (file)
--- a/src/ParaMEDMEM/DEC.cxx
+++ b/src/ParaMEDMEM/DEC.cxx
@@ -188,12 +188,12 @@ namespace ParaMEDMEM
     \f]
 
   */
-  void DEC::renormalizeTargetField()
+  void DEC::renormalizeTargetField(bool isWAbs)
   {
     if (_source_group->containsMyRank())
       for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
         {
-          double total_norm = _local_field->getVolumeIntegral(icomp+1);
+          double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
           double source_norm = total_norm;
           _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 
@@ -202,7 +202,7 @@ namespace ParaMEDMEM
       {
         for (int icomp=0; icomp<_local_field->getField()->getArray()->getNumberOfComponents(); icomp++)
           {
-            double total_norm = _local_field->getVolumeIntegral(icomp+1);
+            double total_norm = _local_field->getVolumeIntegral(icomp+1,isWAbs);
             double source_norm=total_norm;
             _comm_interface->broadcast(&source_norm, 1, MPI_DOUBLE, 0,* dynamic_cast<MPIProcessorGroup*>(_union_group)->getComm());
 

diff --git a/src/ParaMEDMEM/DEC.hxx b/src/ParaMEDMEM/DEC.hxx
index 0f4c545bbe5066d78c83704d02f5804111565853..d6cc17a7bced119b8b6fa030a917c46c61849421 100644 (file)
--- a/src/ParaMEDMEM/DEC.hxx
+++ b/src/ParaMEDMEM/DEC.hxx
@@ -50,7 +50,7 @@ namespace ParaMEDMEM
     virtual void synchronize()=0;
     virtual ~DEC();
     virtual void computeProcGroup() { }
-    void renormalizeTargetField();
+    void renormalizeTargetField(bool isWAbs);
   protected:
     void compareFieldAndMethod() const throw(INTERP_KERNEL::Exception);
   protected:

diff --git a/src/ParaMEDMEM/InterpKernelDEC.cxx b/src/ParaMEDMEM/InterpKernelDEC.cxx
index 6badf51170b51de7d9c29d023cc14d62a937247c..444378a1fbbad76766213e4ecdb92fd707677640 100644 (file)
--- a/src/ParaMEDMEM/InterpKernelDEC.cxx
+++ b/src/ParaMEDMEM/InterpKernelDEC.cxx
@@ -222,7 +222,7 @@ namespace ParaMEDMEM
       {
         _interpolation_matrix->multiply(*_local_field->getField());
         if (getForcedRenormalization())
-          renormalizeTargetField();
+          renormalizeTargetField(getMeasureAbsStatus());
       }
   }
 
@@ -249,7 +249,7 @@ namespace ParaMEDMEM
     
         _interpolation_matrix->multiply(*_local_field->getField());
         if (getForcedRenormalization())
-          renormalizeTargetField();
+          renormalizeTargetField(getMeasureAbsStatus());
     
       }
     else if (_target_group->containsMyRank())

diff --git a/src/ParaMEDMEM/InterpolationMatrix.cxx b/src/ParaMEDMEM/InterpolationMatrix.cxx
index 6edb3d40e1bc35706676133103803481342cf059..42b6e8d635176638c32d9a6dc72401bcf75f96e0 100644 (file)
--- a/src/ParaMEDMEM/InterpolationMatrix.cxx
+++ b/src/ParaMEDMEM/InterpolationMatrix.cxx
@@ -154,7 +154,7 @@ namespace ParaMEDMEM
 
     MEDCouplingFieldDouble *target_triangle_surf=0;
     if(needTargetSurf)
-      target_triangle_surf = distant_support.getMeasureField();
+      target_triangle_surf = distant_support.getMeasureField(getMeasureAbsStatus());
     fillDSFromVM(iproc_distant,distant_elems,surfaces,target_triangle_surf);
 
     if(needTargetSurf)
@@ -320,7 +320,7 @@ namespace ParaMEDMEM
 
   void InterpolationMatrix::computeIntegralDenoW(ElementLocator& elementLocator)
   {
-    MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField();
+    MEDCouplingFieldDouble *source_triangle_surf = _source_support->getMeasureField(getMeasureAbsStatus());
     _deno_multiply.resize(_coeffs.size());
     vector<vector<double> >::iterator iter6=_deno_multiply.begin();
     const double *values=source_triangle_surf->getArray()->getConstPointer();

diff --git a/src/ParaMEDMEM/ParaFIELD.cxx b/src/ParaMEDMEM/ParaFIELD.cxx
index 6afb239116282fc122739a9c113948c66a195cbd..58532760adc6e7be9d0118a5d28e181a34b32a26 100644 (file)
--- a/src/ParaMEDMEM/ParaFIELD.cxx
+++ b/src/ParaMEDMEM/ParaFIELD.cxx
@@ -211,10 +211,10 @@ namespace ParaMEDMEM
 
   /*! This method retrieves the integral of component \a icomp
     over the all domain. */
-  double ParaFIELD::getVolumeIntegral(int icomp) const
+  double ParaFIELD::getVolumeIntegral(int icomp, bool isWAbs) const
   {
     CommInterface comm_interface = _topology->getProcGroup()->getCommInterface();
-    double integral=_field->measureAccumulate(icomp);
+    double integral=_field->measureAccumulate(icomp,isWAbs);
     double total=0.;
     const MPI_Comm* comm = (dynamic_cast<const MPIProcessorGroup*>(_topology->getProcGroup()))->getComm();
     comm_interface.allReduce(&integral, &total, 1, MPI_DOUBLE, MPI_SUM, *comm);

diff --git a/src/ParaMEDMEM/ParaFIELD.hxx b/src/ParaMEDMEM/ParaFIELD.hxx
index a9bcf57b88d4a63af70f9de4f070cef4d7f89c2a..22adbf5c043ca55fda31315129ecf2e900263c94 100644 (file)
--- a/src/ParaMEDMEM/ParaFIELD.hxx
+++ b/src/ParaMEDMEM/ParaFIELD.hxx
@@ -48,7 +48,7 @@ namespace ParaMEDMEM
     Topology* getTopology() const { return _topology; }
     ParaMESH* getSupport() const  { return _support; }
     int nbComponents() const;
-    double getVolumeIntegral(int icomp) const;
+    double getVolumeIntegral(int icomp, bool isWAbs) const;
     double getL2Norm()const { return -1; }
   private:
     MEDCouplingFieldDouble* _field;

diff --git a/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx b/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx
index cfaef00ff69a462886a0e23a08974a6cc2b2850b..ed54396b6ec15edc87bb8ea7a761037d37e433a7 100644 (file)
--- a/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx
+++ b/src/ParaMEDMEM/Test/ParaMEDMEMTest_InterpKernelDEC.cxx
@@ -199,7 +199,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_2D_(const char *srcMeth, const char *ta
   
   if (source_group->containsMyRank())
     { 
-      field_before_int = parafield->getVolumeIntegral(0);
+      field_before_int = parafield->getVolumeIntegral(0,true);
       dec.synchronize();
       cout<<"DEC usage"<<endl;
       dec.setForcedRenormalization(false);
@@ -223,7 +223,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_2D_(const char *srcMeth, const char *ta
      
       filename<<"./sourcesquare_"<<source_group->myRank()+1;
       aRemover.Register(filename.str().c_str());
-      field_after_int = parafield->getVolumeIntegral(0);
+      field_after_int = parafield->getVolumeIntegral(0,true);
       
       
       //      MPI_Bcast(&field_before_int,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
@@ -406,7 +406,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_3D_(const char *srcMeth, const char *ta
   
   if (source_group->containsMyRank())
     { 
-      field_before_int = parafield->getVolumeIntegral(0);
+      field_before_int = parafield->getVolumeIntegral(0,true);
       dec.synchronize();
       cout<<"DEC usage"<<endl;
       dec.setForcedRenormalization(false);
@@ -430,7 +430,7 @@ void ParaMEDMEMTest::testInterpKernelDEC_3D_(const char *srcMeth, const char *ta
      
       filename<<"./sourcesquare_"<<source_group->myRank()+1;
       aRemover.Register(filename.str().c_str());
-      field_after_int = parafield->getVolumeIntegral(0);
+      field_after_int = parafield->getVolumeIntegral(0,true);
 
       CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, 1e-6);
     
@@ -1190,7 +1190,7 @@ void ParaMEDMEMTest::testAsynchronousInterpKernelDEC_2D(double dtA, double tmaxA
           cout << "testAsynchronousInterpKernelDEC_2D" << rank << " time " << time
                << " dtB " << dtB << " tmaxB " << tmaxB << endl ;
           dec.recvData( time );
-          double vi = parafield->getVolumeIntegral(0);
+          double vi = parafield->getVolumeIntegral(0,true);
           cout << "testAsynchronousInterpKernelDEC_2D" << rank << " time " << time
                << " VolumeIntegral " << vi
                << " time*10000 " << time*10000 << endl ;

diff --git a/src/ParaMEDMEM/Test/test_perf.cxx b/src/ParaMEDMEM/Test/test_perf.cxx
index 72ac18b652dc096692e953cd4d76d1bf0489ec1e..f5872c15bb7dc82b48b93275a75afcf73492333b 100644 (file)
--- a/src/ParaMEDMEM/Test/test_perf.cxx
+++ b/src/ParaMEDMEM/Test/test_perf.cxx
@@ -210,7 +210,7 @@ void testInterpKernelDEC_2D(const string& filename_xml1, const string& meshname1
   double field_after_int;
   
   if (source_group->containsMyRank()){ 
-    field_before_int = parafield->getVolumeIntegral(0);
+    field_before_int = parafield->getVolumeIntegral(0,true);
     get_time( &telps, &tcpu_u, &tcpu_s, &tcpu );
     dec.synchronize();
     get_time( &telps, &tcpu_u, &tcpu_s, &tcpu );
@@ -229,7 +229,7 @@ void testInterpKernelDEC_2D(const string& filename_xml1, const string& meshname1
       cout << "SEND DATA : Telapse = " << telps << " TuserCPU = " << tcpu_u << " TsysCPU = " << tcpu_s << " TCPU = " << tcpu << endl;
     dec.recvData();
      
-    field_after_int = parafield->getVolumeIntegral(0);
+    field_after_int = parafield->getVolumeIntegral(0,true);
 //    CPPUNIT_ASSERT_DOUBLES_EQUAL(field_before_int, field_after_int, epsilon);
       
   }